bug 13918: Unmanaged operations on hypermatrix did not call overload functions.
[scilab.git] / scilab / modules / ast / src / cpp / operations / types_divide.cpp
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2008-2008 - DIGITEO - Antoine ELIAS
4  *
5  *  This file must be used under the terms of the CeCILL.
6  *  This source file is licensed as described in the file COPYING, which
7  *  you should have received as part of this distribution.  The terms
8  *  are also available at
9  *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10  *
11  */
12
13 #include "types_divide.hxx"
14 #include "types_finite.hxx"
15
16 extern "C"
17 {
18 #include "matrix_right_division.h"
19 #include "sciprint.h"
20 #include "localization.h"
21 #include "charEncoding.h"
22 #include "configvariable_interface.h"
23 #include "elem_common.h"
24 }
25
26 using namespace types;
27
28 InternalType *GenericRDivide(InternalType *_pLeftOperand, InternalType *_pRightOperand)
29 {
30     InternalType *pResult       = NULL;
31     GenericType::ScilabType TypeL = _pLeftOperand->getType();
32     GenericType::ScilabType TypeR = _pRightOperand->getType();
33
34     int iResult = 0;
35
36     if (_pLeftOperand->isDouble() && _pLeftOperand->getAs<Double>()->isEmpty())
37     {
38         return Double::Empty();
39     }
40
41     if (_pRightOperand->isDouble() && _pRightOperand->getAs<Double>()->isEmpty())
42     {
43         return Double::Empty();
44     }
45
46     /*
47     ** DOUBLE / DOUBLE
48     */
49     if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabDouble)
50     {
51         Double *pL  = _pLeftOperand->getAs<Double>();
52         Double *pR  = _pRightOperand->getAs<Double>();
53
54         iResult = RDivideDoubleByDouble(pL, pR, (Double**)&pResult);
55     }
56
57     /*
58     ** POLY / DOUBLE
59     */
60     else if (TypeL == GenericType::ScilabPolynom && TypeR == GenericType::ScilabDouble)
61     {
62         Polynom *pL = _pLeftOperand->getAs<types::Polynom>();
63         Double *pR  = _pRightOperand->getAs<Double>();
64
65         iResult = RDividePolyByDouble(pL, pR, (Polynom**)&pResult);
66     }
67
68     /*
69     ** DOUBLE / POLY
70     */
71     else if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabPolynom)
72     {
73         Double *pL  = _pLeftOperand->getAs<Double>();
74         Polynom *pR = _pRightOperand->getAs<types::Polynom>();
75
76         iResult = RDivideDoubleByPoly(pL, pR, (Polynom**)&pResult);
77     }
78
79     /*
80     ** SPARSE / DOUBLE
81     */
82     else if (TypeL == GenericType::ScilabSparse && TypeR == GenericType::ScilabDouble)
83     {
84         Sparse *pL = _pLeftOperand->getAs<types::Sparse>();
85         Double *pR = _pRightOperand->getAs<Double>();
86
87         iResult = RDivideSparseByDouble(pL, pR, &pResult);
88     }
89
90     //manage errors
91     if (iResult)
92     {
93         switch (iResult)
94         {
95             case 1 :
96                 throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
97             case 2 :
98                 throw ast::InternalError(_W("With NaN or Inf a division by scalar expected.\n"));
99             case 3 :
100                 throw ast::InternalError(_W("Division by zero...\n"));
101             case 4 :
102                 if (getWarningMode())
103                 {
104                     sciprint(_("Warning : Division by zero...\n"));
105                 }
106                 break;
107             //            default : throw ast::InternalError(_W("Operator / : Error %d not yet managed.\n"), iResult);
108             default :
109                 sciprint(_("Operator / : Error %d not yet managed.\n"), iResult);
110         }
111     }
112
113     /*
114     ** Default case : Return NULL will Call Overloading.
115     */
116     return pResult;
117 }
118
119
120 int RDivideDoubleByDouble(Double *_pDouble1, Double *_pDouble2, Double **_pDoubleOut)
121 {
122     int iErr = 0;
123
124     //check finite values of _pDouble1 and _pDouble2
125     if (isDoubleFinite(_pDouble1) == false || isDoubleFinite(_pDouble2) == false)
126     {
127         if (_pDouble2->isScalar() == false)
128         {
129             return 2;
130         }
131     }
132
133     if (_pDouble2->isScalar())
134     {
135         //Y / x
136         int iInc1       = 1;
137         int iInc2       = 0;
138         bool bComplex1  = _pDouble1->isComplex();
139         bool bComplex2  = _pDouble2->isComplex();
140
141         *_pDoubleOut    = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2);
142
143         if (bComplex1 == false && bComplex2 == false)
144         {
145             // Real1 \ Real2 -> Real2 / Real1
146             iErr = iRightDivisionRealMatrixByRealMatrix(
147                        _pDouble1->get(), iInc1,
148                        _pDouble2->get(), iInc2,
149                        (*_pDoubleOut)->get(), 1, _pDouble1->getSize());
150         }
151         else if (bComplex1 == false && bComplex2 == true)
152         {
153             // Real \ Complex -> Complex / Real
154             iErr = iRightDivisionRealMatrixByComplexMatrix(
155                        _pDouble1->get(), iInc1,
156                        _pDouble2->get(), _pDouble2->getImg(), iInc2,
157                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
158         }
159         else if (bComplex1 == true && bComplex2 == false)
160         {
161             // Complex \ Real -> Real / Complex
162             iErr = iRightDivisionComplexMatrixByRealMatrix(
163                        _pDouble1->get(), _pDouble1->getImg(), iInc1,
164                        _pDouble2->get(), iInc2,
165                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
166         }
167         else if (bComplex1 == true && bComplex2 == true)
168         {
169             // Complex \ Complex
170             iErr = iRightDivisionComplexMatrixByComplexMatrix(
171                        _pDouble1->get(), _pDouble1->getImg(), iInc1,
172                        _pDouble2->get(), _pDouble2->getImg(), iInc2,
173                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
174         }
175
176         return iErr;
177     }
178
179     if (_pDouble1->isScalar())
180     {
181         if (_pDouble2->getDims() > 2)
182         {
183             //not managed, call overload
184             return 0;
185         }
186
187         // x / eye() = x
188         if (_pDouble2->isIdentity() )
189         {
190             *_pDoubleOut    = new Double(*_pDouble1);
191             return 0;
192         }
193         double dblSavedR = 0;
194         double dblSavedI = 0;
195         Double *pdblTemp = NULL;
196
197         int iRowResult = _pDouble2->getCols();
198         int iColResult = _pDouble2->getRows();
199
200
201
202
203         //in this case, we have to create a temporary square matrix
204         pdblTemp = new Double(iRowResult, iRowResult, _pDouble1->isComplex());
205         pdblTemp->setZeros();
206
207         if (_pDouble1->isComplex())
208         {
209             dblSavedR = _pDouble1->getReal()[0];
210             dblSavedI = _pDouble1->getImg()[0];
211             for (int i = 0 ; i < iRowResult ; i++)
212             {
213                 pdblTemp->set(i, i, dblSavedR);
214                 pdblTemp->setImg(i, i, dblSavedI);
215             }
216         }
217         else
218         {
219             dblSavedR = _pDouble1->getReal()[0];
220             for (int i = 0 ; i < iRowResult ; i++)
221             {
222                 pdblTemp->set(i, i, dblSavedR);
223             }
224         }
225
226         *_pDoubleOut = new Double(iRowResult, iColResult, _pDouble1->isComplex() || _pDouble2->isComplex());
227
228         if ((*_pDoubleOut)->isComplex())
229         {
230             double dblRcond = 0;
231             iErr = iRightDivisionOfComplexMatrix(
232                        pdblTemp->getReal(), pdblTemp->getImg(), pdblTemp->getRows(), pdblTemp->getCols(),
233                        _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
234                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iRowResult, iColResult, &dblRcond);
235         }
236         else
237         {
238             double dblRcond = 0;
239             iErr = iRightDivisionOfRealMatrix(
240                        pdblTemp->getReal(), pdblTemp->getRows(), pdblTemp->getCols(),
241                        _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
242                        (*_pDoubleOut)->getReal(), iRowResult, iColResult, &dblRcond);
243         }
244         delete pdblTemp;
245         return iErr;
246     }
247
248     if (_pDouble1->getDims() > 2 || _pDouble2->getDims() > 2 || _pDouble1->getCols() != _pDouble2->getCols())
249     {
250         //not managed
251         return 0;
252     }
253
254     *_pDoubleOut = new Double(_pDouble1->getRows(), _pDouble2->getRows(), _pDouble1->isComplex() || _pDouble2->isComplex());
255     if ((*_pDoubleOut)->isComplex())
256     {
257         double dblRcond = 0;
258         iErr = iRightDivisionOfComplexMatrix(
259                    _pDouble1->getReal(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
260                    _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
261                    (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), _pDouble1->getRows(), _pDouble2->getRows(), &dblRcond);
262     }
263     else
264     {
265         double dblRcond = 0;
266         iErr = iRightDivisionOfRealMatrix(
267                    _pDouble1->getReal(), _pDouble1->getRows(), _pDouble1->getCols(),
268                    _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
269                    (*_pDoubleOut)->getReal(), _pDouble1->getRows(), _pDouble2->getRows(), &dblRcond);
270     }
271
272     return iErr;
273 }
274
275 int RDividePolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom** _pPolyOut)
276 {
277     bool bComplex1  = _pPoly->isComplex();
278     bool bComplex2  = _pDouble->isComplex();
279     bool bScalar1   = _pPoly->getRows() == 1  && _pPoly->getCols() == 1;
280     bool bScalar2   = _pDouble->getRows() == 1 && _pDouble->getCols() == 1;
281
282     Polynom *pTemp = NULL; //use only if _pPoly is scalar and _pDouble not.
283
284     int iRowResult  = 0;
285     int iColResult = 0;
286     int *piRank   = NULL;
287
288     /* if(bScalar1 && bScalar2)
289     {
290     iRowResult = 1;
291     iColResult = 1;
292
293     piRank = new int[1];
294     piRank[0] = _pPoly->get(0)->getRank();
295     }
296     else */
297
298     if (bScalar1 == false && bScalar2 == false)
299     {
300         // call overload
301         return 0;
302     }
303
304     if (bScalar2)
305     {
306         double dblDivR = _pDouble->get(0);
307         double dblDivI = _pDouble->getImg(0);
308
309         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
310         if (_pDouble->isComplex())
311         {
312             (*_pPolyOut)->setComplex(true);
313         }
314
315         for (int i = 0 ; i < _pPoly->getSize() ; i++)
316         {
317             bool bComplex1 = _pPoly->isComplex();
318             bool bComplex2 = _pDouble->isComplex();
319
320             SinglePoly* pC = (*_pPolyOut)->get(i);
321
322             if (bComplex1 == false && bComplex2 == false)
323             {
324                 iRightDivisionRealMatrixByRealMatrix(pC->get(), 1, &dblDivR, 0, pC->get(), 1, pC->getSize());
325             }
326             else if (bComplex1 == true && bComplex2 == false)
327             {
328                 iRightDivisionComplexMatrixByRealMatrix(pC->get(), pC->getImg(), 1, &dblDivR, 0, pC->get(), pC->getImg(), 1, pC->getSize());
329             }
330             else if (bComplex1 == false && bComplex2 == true)
331             {
332                 iRightDivisionRealMatrixByComplexMatrix(pC->get(), 1, &dblDivR, &dblDivI, 0, pC->get(), pC->getImg(), 1, pC->getSize());
333             }
334             else if (bComplex1 == true && bComplex2 == true)
335             {
336                 iRightDivisionComplexMatrixByComplexMatrix(pC->get(), pC->getImg(), 1, &dblDivR, &dblDivI, 0, pC->get(), pC->getImg(), 1, pC->getSize());
337             }
338         }
339
340         return 0;
341     }
342
343     if (bScalar1)
344     {
345         //in this case, we have to create a temporary square polinomial matrix
346         iRowResult = _pDouble->getCols();
347         iColResult = _pDouble->getRows();
348
349         piRank = new int[iRowResult * iRowResult];
350         int iMaxRank = _pPoly->getMaxRank();
351         for (int i = 0 ; i < iRowResult * iRowResult ; i++)
352         {
353             piRank[i] = iMaxRank;
354         }
355
356         pTemp = new Polynom(_pPoly->getVariableName(), iRowResult, iRowResult, piRank);
357         if (bComplex1 || bComplex2)
358         {
359             pTemp->setComplex(true);
360         }
361
362         SinglePoly *pdblData = _pPoly->get(0);
363         for (int i = 0 ; i < iRowResult ; i++)
364         {
365             pTemp->set(i, i, pdblData);
366         }
367     }
368
369     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), iRowResult, iColResult, piRank);
370     delete[] piRank;
371     if (bComplex1 || bComplex2)
372     {
373         (*_pPolyOut)->setComplex(true);
374     }
375
376     if (bScalar2)
377     {
378         //[p] * cst
379         for (int i = 0 ; i < _pPoly->getSize() ; i++)
380         {
381             SinglePoly *pPolyIn   = _pPoly->get(i);
382             double* pRealIn  = pPolyIn->get();
383             double* pImgIn  = pPolyIn->getImg();
384
385             SinglePoly *pPolyOut  = (*_pPolyOut)->get(i);
386             double* pRealOut = pPolyOut->get();
387             double* pImgOut  = pPolyOut->getImg();
388
389             if (bComplex1 == false && bComplex2 == false)
390             {
391                 iRightDivisionRealMatrixByRealMatrix(
392                     pRealIn, 1,
393                     _pDouble->getReal(), 0,
394                     pRealOut, 1, pPolyOut->getSize());
395             }
396             else if (bComplex1 == false && bComplex2 == true)
397             {
398                 iRightDivisionRealMatrixByComplexMatrix(
399                     pRealIn, 1,
400                     _pDouble->getReal(), _pDouble->getImg(), 0,
401                     pRealOut, pImgOut, 1, pPolyOut->getSize());
402             }
403             else if (bComplex1 == true && bComplex2 == false)
404             {
405                 iRightDivisionComplexMatrixByRealMatrix(
406                     pRealIn, pImgIn, 1,
407                     _pDouble->getReal(), 0,
408                     pRealOut, pImgOut, 1, pPolyOut->getSize());
409             }
410             else if (bComplex1 == true && bComplex2 == true)
411             {
412                 iRightDivisionComplexMatrixByComplexMatrix(
413                     pRealIn, pImgIn, 1,
414                     _pDouble->getReal(), _pDouble->getImg(), 0,
415                     pRealOut, pImgOut, 1, pPolyOut->getSize());
416             }
417         }
418     }
419     else if (bScalar1)
420     {
421         for (int i = 0 ; i < pTemp->get(0)->getSize() ; i++)
422         {
423             Double *pCoef    = pTemp->extractCoef(i);
424             Double *pResultCoef = new Double(iRowResult, iColResult, pCoef->isComplex());
425             double *pReal    = pResultCoef->getReal();
426             double *pImg    = pResultCoef->getImg();
427
428             if (bComplex1 == false && bComplex2 == false)
429             {
430                 double dblRcond = 0;
431                 iRightDivisionOfRealMatrix(
432                     pCoef->getReal(), iRowResult, iRowResult,
433                     _pDouble->getReal(), _pDouble->getRows(), _pDouble->getCols(),
434                     pReal, iRowResult, iColResult, &dblRcond);
435             }
436             else
437             {
438                 double dblRcond = 0;
439                 iRightDivisionOfComplexMatrix(
440                     pCoef->getReal(), pCoef->getImg(), iRowResult, iRowResult,
441                     _pDouble->getReal(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
442                     pReal, pImg, iRowResult, iColResult, &dblRcond);
443             }
444
445             (*_pPolyOut)->insertCoef(i, pResultCoef);
446             delete pCoef;
447             delete pResultCoef;
448         }
449     }
450     return 0;
451 }
452
453 int RDivideDoubleByPoly(Double* /*_pDouble*/, Polynom* /*_pPoly*/, Polynom** /*_pPolyOut*/)
454 {
455     return 0;
456 }
457
458 int RDivideSparseByDouble(types::Sparse* _pSp, types::Double* _pDouble, InternalType** _pSpOut)
459 {
460     if (_pDouble->isEmpty())
461     {
462         //sp / []
463         *_pSpOut = Double::Empty();
464         return 0;
465     }
466
467     if (_pDouble->isIdentity())
468     {
469         *_pSpOut    = new Sparse(*_pSp);
470         return 0;
471     }
472
473     size_t iSize = _pSp->nonZeros();
474     int* Col = new int[iSize];
475     int* Row = new int[_pSp->getRows()];
476     _pSp->getColPos(Col);
477     _pSp->getNbItemByRow(Row);
478     int* iPositVal = new int[iSize];
479
480     int idx = 0;
481     for (int i = 0; i < _pSp->getRows(); i++)
482     {
483         for (int j = 0; j < Row[i]; j++)
484         {
485             iPositVal[idx] = (Col[idx] - 1) * _pSp->getRows() + i;
486             ++idx;
487         }
488     }
489
490     delete[] Col;
491     delete[] Row;
492
493     Double** pDbl = new Double*[iSize];
494     Double** pDblSp = new Double*[iSize];
495     double* pdbl = _pDouble->get();
496
497     if (_pDouble->isScalar())
498     {
499         if (_pDouble->isComplex())
500         {
501             double* pdblImg = _pDouble->getImg();
502             for (int i = 0; i < iSize; i++)
503             {
504                 pDbl[i] = new Double(pdbl[0], pdblImg[0]);
505                 pDblSp[i] = new Double(_pSp->get(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
506             }
507         }
508         else
509         {
510             for (int i = 0; i < iSize; i++)
511             {
512                 pDbl[i] = new Double(pdbl[0]);
513                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
514             }
515         }
516     }
517     else if (_pDouble->getSize() == iSize)
518     {
519         if (_pDouble->isComplex())
520         {
521             double* pdblImg = _pDouble->getImg();
522             for (int i = 0; i < iSize; i++)
523             {
524                 pDbl[i] = new Double(pdbl[i], pdblImg[i]);
525                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
526             }
527         }
528         else
529         {
530             for (int i = 0; i < iSize; i++)
531             {
532                 pDbl[i] = new Double(pdbl[i]);
533                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
534             }
535         }
536     }
537     else
538     {
539         delete[] pDbl;
540         delete[] pDblSp;
541         delete[] iPositVal;
542         return 0;
543     }
544
545     Sparse* pSpTemp = new Sparse(_pSp->getRows(), _pSp->getCols(), _pSp->isComplex() || _pDouble->isComplex());
546     pSpTemp->zero_set();
547
548     Double* ppDblGet = NULL;
549     int iResultat;
550     for (int i = 0; i < iSize; i++)
551     {
552         if ((pDblSp[i]->get(0) != 0) || (pDblSp[i]->getImg(0) != 0))
553         {
554             iResultat = RDivideDoubleByDouble(pDblSp[i], pDbl[i], &ppDblGet);
555             if (iResultat != 0)
556             {
557                 delete ppDblGet;
558                 return iResultat;
559             }
560             std::complex<double> cplx(ppDblGet->get(0), ppDblGet->getImg(0));
561             pSpTemp->set(iPositVal[i], cplx, true);
562             delete ppDblGet;
563         }
564     }
565
566     delete[] iPositVal;
567
568     for (int i = 0; i < iSize; ++i)
569     {
570         delete pDbl[i];
571         delete pDblSp[i];
572     }
573
574     delete[] pDbl;
575     delete[] pDblSp;
576
577     *_pSpOut = pSpTemp;
578     return 0;
579 }