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