Coverity #1247398, #1199692, #1321278, #1321279, #1321280, #1321281, #1097625 fixed...
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_spec.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2009 - DIGITEO - Bernard HUGUENEY
4 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
5 *
6  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14 *
15 */
16 /*--------------------------------------------------------------------------*/
17
18 #include "linear_algebra_gw.hxx"
19 #include "function.hxx"
20 #include "double.hxx"
21 #include "overload.hxx"
22
23 extern "C"
24 {
25 #include "localization.h"
26 #include "Scierror.h"
27 #include "sciprint.h"
28 #include "eigen.h"
29 #include "issymmetric.h"
30 #include "vfinite.h"
31 }
32
33 bool isNoZeroImag(types::Double* _pDbl);
34 /*--------------------------------------------------------------------------*/
35
36 types::Function::ReturnValue sci_spec(types::typed_list &in, int _iRetCount, types::typed_list &out)
37 {
38     double* pDataA          = NULL;
39     double* pDataB          = NULL;
40     bool symmetric          = FALSE;
41     int iRet                = 0;
42
43     if (in.size() != 1 && in.size() != 2)
44     {
45         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "spec", 1, 2);
46         return types::Function::Error;
47     }
48
49     if (_iRetCount > 2 * in.size())
50     {
51         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "spec", 1, 2 * in.size());
52         return types::Function::Error;
53     }
54
55     if (in[0]->isDouble() == false)
56     {
57         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_spec";
58         return Overload::call(wstFuncName, in, _iRetCount, out);
59     }
60
61     types::Double* in0 = in[0]->getAs<types::Double>();
62
63     if (in0->getCols() != in0->getRows())
64     {
65         Scierror(20, _("%s: Wrong type for input argument #%d: Square matrix expected.\n"), "spec", 1);
66         return types::Function::Error;
67     }
68
69     if (in0->getRows() == -1 || in0->getCols() == -1) // manage eye case
70     {
71         Scierror(271, _("%s: Size varying argument a*eye(), (arg %d) not allowed here.\n"), "spec", 1);
72         return types::Function::Error;
73     }
74
75     if (in0->getCols() == 0 || in0->getRows() == 0) // size null
76     {
77         out.push_back(types::Double::Empty());
78         for (int i = 1; i < _iRetCount; i++)
79         {
80             out.push_back(types::Double::Empty());
81         }
82         return types::Function::OK;
83     }
84
85     types::Double* pDblA = in0->clone()->getAs<types::Double>();
86
87     if (in.size() == 1)
88     {
89         types::Double* pDblEigenValues  = NULL;
90         types::Double* pDblEigenVectors = NULL;
91
92         if (pDblA->isComplex())
93         {
94             pDataA = (double*)oGetDoubleComplexFromPointer(pDblA->getReal(), pDblA->getImg(), pDblA->getSize());
95             if (!pDataA)
96             {
97                 pDblA->killMe();
98                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
99                 return types::Function::Error;
100             }
101         }
102         else
103         {
104             pDataA = pDblA->getReal();
105         }
106
107         int totalSize = pDblA->getSize();
108         if ((pDblA->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataA) : C2F(vfinite)(&totalSize, pDataA)) == false)
109         {
110             if (pDblA->isComplex())
111             {
112                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
113             }
114             pDblA->killMe();
115             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 1);
116             return types::Function::Error;
117         }
118
119         symmetric = isSymmetric(pDblA->getReal(), pDblA->getImg(), pDblA->isComplex(), pDblA->getRows(), pDblA->getCols()) == 1;
120         int eigenValuesCols = (_iRetCount == 1) ? 1 : pDblA->getCols();
121
122         if (symmetric)
123         {
124             pDblEigenValues = new types::Double(pDblA->getCols(), eigenValuesCols);
125             if (_iRetCount == 2)
126             {
127                 pDblEigenVectors = new types::Double(pDblA->getCols(), pDblA->getCols(), pDblA->isComplex());
128             }
129         }
130         else
131         {
132             pDblEigenValues  = new types::Double(pDblA->getCols(), eigenValuesCols, true);
133             if (_iRetCount == 2)
134             {
135                 pDblEigenVectors = new types::Double(pDblA->getCols(), pDblA->getCols(), true);
136             }
137         }
138
139         if (pDblA->isComplex())
140         {
141             if (symmetric)
142             {
143                 iRet = iEigen1ComplexSymmetricM((doublecomplex*)pDataA, pDblA->getCols(), (_iRetCount == 2), pDblEigenValues->getReal());
144
145                 if (iRet < 0)
146                 {
147                     vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
148                     pDblA->killMe();
149                     Scierror(998, _("%s: On entry to ZGEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
150                     return types::Function::Error;
151                 }
152
153                 if (iRet > 0)
154                 {
155                     vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
156                     pDblA->killMe();
157                     Scierror(24, _("%s: Convergence problem, %d off-diagonal elements of an intermediate tridiagonal form did not converge to zero.\n"), "spec", iRet);
158                     return types::Function::Error;
159                 }
160
161                 if (_iRetCount == 2)
162                 {
163                     vGetPointerFromDoubleComplex((doublecomplex*)pDataA, pDblA->getSize() , pDblEigenVectors->getReal(), pDblEigenVectors->getImg());
164                     vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
165                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
166                     out.push_back(pDblEigenVectors);
167                 }
168                 out.push_back(pDblEigenValues);
169                 pDblA->killMe();
170             }
171             else // not symmetric
172             {
173                 doublecomplex* pEigenValues = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
174                 doublecomplex* pEigenVectors = pDblEigenVectors ? (doublecomplex*)MALLOC(sizeof(doublecomplex) * pDblA->getSize()) : NULL;
175                 iRet = iEigen1ComplexM((doublecomplex*)pDataA, pDblA->getCols(), pEigenValues, pEigenVectors);
176                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
177                 if (iRet < 0)
178                 {
179                     pDblA->killMe();
180                     Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
181                     return types::Function::Error;
182                 }
183
184                 if (iRet > 0)
185                 {
186                     pDblA->killMe();
187                     Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of W contain eigenvalues which have converged.\n"), "spec", iRet);
188                     return types::Function::Error;
189                 }
190
191                 if (_iRetCount == 2)
192                 {
193                     expandZToDiagonalOfCMatrix(pEigenValues, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg());
194                     vGetPointerFromDoubleComplex(pEigenVectors, pDblA->getSize(), pDblEigenVectors->getReal(), pDblEigenVectors->getImg());
195
196                     FREE(pEigenVectors);
197                     out.push_back(pDblEigenVectors);
198                 }
199                 else
200                 {
201                     vGetPointerFromDoubleComplex(pEigenValues, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg());
202                 }
203                 out.push_back(pDblEigenValues);
204                 FREE(pEigenValues);
205                 pDblA->killMe();
206             }
207         }
208         else // real
209         {
210             if (symmetric)
211             {
212                 iRet = iEigen1RealSymmetricM(pDataA, pDblA->getCols(), (_iRetCount == 2), pDblEigenValues->getReal());
213                 if (iRet < 0)
214                 {
215                     pDblA->killMe();
216                     Scierror(998, _("%s: On entry to ZGEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
217                     return types::Function::Error;
218                 }
219
220                 if (iRet > 0)
221                 {
222                     pDblA->killMe();
223                     Scierror(24, _("%s: Convergence problem, %d off-diagonal elements of an intermediate tridiagonal form did not converge to zero.\n"), "spec", iRet);
224                     return types::Function::Error;
225                 }
226
227                 if (_iRetCount == 2)
228                 {
229                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
230                     out.push_back(pDblA);
231                 }
232                 else
233                 {
234                     pDblA->killMe();
235                 }
236
237                 out.push_back(pDblEigenValues);
238             }
239             else // not symmetric
240             {
241                 iRet = iEigen1RealM(pDataA, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg(), pDblEigenVectors ? pDblEigenVectors->getReal() : NULL, pDblEigenVectors ? pDblEigenVectors->getImg() : NULL);
242
243                 if (iRet < 0)
244                 {
245                     pDblA->killMe();
246                     Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
247                     return types::Function::Error;
248                 }
249
250                 if (iRet > 0)
251                 {
252                     pDblA->killMe();
253                     Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of WR and WI contain eigenvalues which have converged.\n"), "spec", iRet);
254                     return types::Function::Error;
255                 }
256
257                 if (_iRetCount == 2)
258                 {
259                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
260                     expandToDiagonalOfMatrix(pDblEigenValues->getImg(), pDblA->getCols());
261                     out.push_back(pDblEigenVectors);
262                 }
263
264                 out.push_back(pDblEigenValues);
265                 pDblA->killMe();
266             }
267         }
268
269         return types::Function::OK;
270     }
271
272     if (in.size() == 2)
273     {
274         types::Double* pDblL            = NULL;
275         types::Double* pDblR            = NULL;
276         types::Double* pDblBeta         = NULL;
277         types::Double* pDblAlpha        = NULL;
278         doublecomplex* pL               = NULL;
279         doublecomplex* pR               = NULL;
280         doublecomplex* pBeta            = NULL;
281         doublecomplex* pAlpha           = NULL;
282         bool bIsComplex                 = false;
283
284         if (in[1]->isDouble() == false)
285         {
286             std::wstring wstFuncName = L"%" + in[1]->getShortTypeStr() + L"_spec";
287             return Overload::call(wstFuncName, in, _iRetCount, out);
288         }
289
290         types::Double* in1 = in[1]->getAs<types::Double>();
291
292         if (in1->getCols() != in1->getRows())
293         {
294             Scierror(20, _("%s: Wrong type for input argument #%d: Square matrix expected.\n"), "spec", 2);
295             return types::Function::Error;
296         }
297
298         if (pDblA->getRows() != in1->getRows() && pDblA->getCols() != in1->getCols())
299         {
300             pDblA->killMe();
301             Scierror(999, _("%s: Arguments %d and %d must have equal dimensions.\n"), "spec", 1, 2);
302             return types::Function::Error;
303         }
304
305         //chekc if A and B are real complex or with imag part at 0
306         if (isNoZeroImag(pDblA) == false && isNoZeroImag(in1) == false)
307         {
308             //view A and B as real matrix
309             bIsComplex = false;
310         }
311         else
312         {
313             bIsComplex = pDblA->isComplex() || in1->isComplex();
314         }
315
316         types::Double* pDblB = in1->clone()->getAs<types::Double>();
317         if (bIsComplex)
318         {
319             if (pDblA->isComplex() == false)
320             {
321                 pDblA->setComplex(true);
322             }
323
324             if (pDblB->isComplex() == false)
325             {
326                 pDblB->setComplex(true);
327             }
328
329             pDataA = (double*)oGetDoubleComplexFromPointer(pDblA->getReal(), pDblA->getImg(), pDblA->getSize());
330             pDataB = (double*)oGetDoubleComplexFromPointer(pDblB->getReal(), pDblB->getImg(), pDblB->getSize());
331
332             if (!pDataA && !pDataB)
333             {
334                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
335                 return types::Function::Error;
336             }
337
338             if (!pDataA)
339             {
340                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataB);
341                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
342                 return types::Function::Error;
343             }
344             
345             if (!pDataB)
346             {
347                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
348                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
349                 return types::Function::Error;
350             }
351         }
352         else
353         {
354             pDataA = pDblA->getReal();
355             pDataB = pDblB->getReal();
356         }
357
358         int totalSize = pDblA->getSize();
359
360         if ((pDblA->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataA) : C2F(vfinite)(&totalSize, pDataA)) == false)
361         {
362             pDblA->killMe();
363             pDblB->killMe();
364             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 1);
365             return types::Function::Error;
366         }
367
368         if ((pDblB->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataB) : C2F(vfinite)(&totalSize, pDataB)) == false)
369         {
370             pDblA->killMe();
371             pDblB->killMe();
372             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 2);
373             return types::Function::Error;
374         }
375
376         switch (_iRetCount)
377         {
378             case 4:
379             {
380                 pDblL = new types::Double(pDblA->getRows(), pDblA->getCols(), true);
381                 if (bIsComplex)
382                 {
383                     pL = (doublecomplex*)MALLOC(pDblA->getSize() * sizeof(doublecomplex));
384                 }
385             }
386             case 3:
387             {
388                 pDblR = new types::Double(pDblA->getRows(), pDblA->getCols(), true);
389                 if (bIsComplex)
390                 {
391                     pR = (doublecomplex*)MALLOC(pDblA->getSize() * sizeof(doublecomplex));
392                 }
393             }
394             case 2:
395             {
396                 if (bIsComplex)
397                 {
398                     pBeta = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
399                 }
400                 pDblBeta = new types::Double(pDblA->getCols(), 1, pBeta ? true : false);
401             }
402             default : // case 1:
403             {
404                 if (bIsComplex)
405                 {
406                     pAlpha = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
407                 }
408                 pDblAlpha = new types::Double(pDblA->getCols(), 1, true);
409             }
410         }
411
412         if (bIsComplex)
413         {
414             iRet = iEigen2ComplexM((doublecomplex*)pDataA, (doublecomplex*)pDataB, pDblA->getCols(), pAlpha, pBeta, pR, pL);
415         }
416         else
417         {
418             iRet = iEigen2RealM(    pDataA, pDataB, pDblA->getCols(),
419                                     pDblAlpha->getReal(), pDblAlpha->getImg(),
420                                     pDblBeta ? pDblBeta->getReal()  : NULL,
421                                     pDblR    ? pDblR->getReal()     : NULL,
422                                     pDblR    ? pDblR->getImg()      : NULL,
423                                     pDblL    ? pDblL->getReal()     : NULL,
424                                     pDblL    ? pDblL->getImg()      : NULL);
425         }
426
427         if (iRet > 0)
428         {
429             sciprint(_("Warning :\n"));
430             sciprint(_("Non convergence in the QZ algorithm.\n"));
431             sciprint(_("The top %d  x %d blocks may not be in generalized Schur form.\n"), iRet);
432         }
433
434         if (iRet < 0)
435         {
436             pDblA->killMe();
437             pDblB->killMe();
438             Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
439             return types::Function::Error;
440         }
441
442         if (iRet > 0)
443         {
444             if (bIsComplex)
445             {
446                 if (iRet <= pDblA->getCols())
447                 {
448                     Scierror(24, _("%s: The QZ iteration failed in DGGEV.\n"), "spec");
449                 }
450                 else
451                 {
452                     if (iRet == pDblA->getCols() + 1)
453                     {
454                         Scierror(999, _("%s: Other than QZ iteration failed in DHGEQZ.\n"), "spec");
455                     }
456                     if (iRet == pDblA->getCols() + 2)
457                     {
458                         Scierror(999, _("%s: Error return from DTGEVC.\n"), "spec");
459                     }
460                 }
461             }
462             else
463             {
464                 Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of W contain eigenvalues which have converged.\n"), "spec", iRet);
465             }
466
467             pDblA->killMe();
468             pDblB->killMe();
469             if (pDataA)
470             {
471                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
472             }
473
474             if (pDataB)
475             {
476                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataB);
477             }
478             return types::Function::Error;
479         }
480
481         if (bIsComplex)
482         {
483             switch (_iRetCount)
484             {
485                 case 4:
486                     vGetPointerFromDoubleComplex(pL, pDblA->getSize(), pDblL->getReal(), pDblL->getImg());
487                 case 3:
488                     vGetPointerFromDoubleComplex(pR, pDblA->getSize(), pDblR->getReal(), pDblR->getImg());
489                 case 2:
490                     vGetPointerFromDoubleComplex(pBeta, pDblA->getCols(), pDblBeta->getReal(), pDblBeta->getImg());
491                 default : // case 1:
492                     vGetPointerFromDoubleComplex(pAlpha, pDblA->getCols(), pDblAlpha->getReal(), pDblAlpha->getImg());
493             }
494         }
495
496         switch (_iRetCount)
497         {
498             case 1:
499             {
500                 out.push_back(pDblAlpha);
501                 break;
502             }
503             case 2:
504             {
505                 out.push_back(pDblAlpha);
506                 out.push_back(pDblBeta);
507                 break;
508             }
509             case 3:
510             {
511                 out.push_back(pDblAlpha);
512                 out.push_back(pDblBeta);
513                 out.push_back(pDblR);
514                 break;
515             }
516             case 4:
517             {
518                 out.push_back(pDblAlpha);
519                 out.push_back(pDblBeta);
520                 out.push_back(pDblL);
521                 out.push_back(pDblR);
522             }
523         }
524
525         if (pAlpha)
526         {
527             vFreeDoubleComplexFromPointer(pAlpha);
528         }
529         if (pBeta)
530         {
531             vFreeDoubleComplexFromPointer(pBeta);
532         }
533         if (pL)
534         {
535             vFreeDoubleComplexFromPointer(pL);
536         }
537         if (pR)
538         {
539             vFreeDoubleComplexFromPointer(pR);
540         }
541         if (bIsComplex && pDblB->isComplex())
542         {
543             vFreeDoubleComplexFromPointer((doublecomplex*)pDataB);
544         }
545         pDblB->killMe();
546
547     } // if(in.size() == 2)
548
549     if (pDblA->isComplex())
550     {
551         vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
552     }
553
554     return types::Function::OK;
555 }
556 /*--------------------------------------------------------------------------*/
557 bool isNoZeroImag(types::Double* _pDbl)
558 {
559     double* pdbl = _pDbl->getImg();
560     if (pdbl)
561     {
562         for (int i = 0 ; i < _pDbl->getSize() ; i++)
563         {
564             if (pdbl[i])
565             {
566                 return true;
567             }
568         }
569     }
570     return false;
571 }
572 /*--------------------------------------------------------------------------*/