* Bug 15321: lu() was leaking memory
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_lu.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
4 *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13 *
14 */
15 /*--------------------------------------------------------------------------*/
16
17 #include "linear_algebra_gw.hxx"
18 #include "function.hxx"
19 #include "double.hxx"
20 #include "overload.hxx"
21
22 extern "C"
23 {
24 #include "localization.h"
25 #include "Scierror.h"
26 #include "lu.h"
27 #include "doublecomplex.h"
28 }
29 /*--------------------------------------------------------------------------*/
30
31 types::Function::ReturnValue sci_lu(types::typed_list &in, int _iRetCount, types::typed_list &out)
32 {
33     types::Double* pDbl     = NULL;
34     types::Double* pDblL    = NULL;
35     types::Double* pDblU    = NULL;
36     types::Double* pDblE    = NULL;
37     double* pdL             = NULL;
38     double* pdU             = NULL;
39     double* pData           = NULL;
40     int iMinRowsCols        = 0;
41
42     if (in.size() != 1)
43     {
44         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "lu", 1);
45         return types::Function::Error;
46     }
47
48     if (_iRetCount < 2 || _iRetCount > 3)
49     {
50         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "lu", 2, 3);
51         return types::Function::Error;
52     }
53
54     if ((in[0]->isDouble() == false))
55     {
56         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lu";
57         return Overload::call(wstFuncName, in, _iRetCount, out);
58     }
59
60     pDbl = in[0]->getAs<types::Double>()->clone()->getAs<types::Double>();
61
62     if (pDbl->isComplex())
63     {
64         pData = (double*)oGetDoubleComplexFromPointer(pDbl->getReal(), pDbl->getImg(), pDbl->getSize());
65         if (!pData)
66         {
67             Scierror(999, _("%s: Cannot allocate more memory.\n"), "lu");
68             pDbl->killMe();
69             return types::Function::Error;
70         }
71     }
72     else
73     {
74         pData = pDbl->getReal();
75     }
76
77     if ((pDbl->getCols() == 0) || (pDbl->getRows() == 0))
78     {
79         out.push_back(types::Double::Empty());
80         out.push_back(types::Double::Empty());
81         if (_iRetCount == 3)
82         {
83             out.push_back(types::Double::Empty());
84         }
85         pDbl->killMe();
86         return types::Function::OK;
87     }
88
89     if ((pDbl->getRows() == -1) || (pDbl->getCols() == -1)) // Rhs(1)=k*eye() => Lhs(1)=eye()
90     {
91         // Lhs(2)=k*eye(), Lhs(3)=eye()
92         pDblL = new types::Double(-1, -1, pDbl->isComplex());
93         pDblL->set(0, 1);
94         out.push_back(pDblL);
95         out.push_back(pDbl);
96         if (_iRetCount == 3)
97         {
98             pDblE = new types::Double(-1, -1, pDbl->isComplex());
99             pDblE->set(0, 1);
100             out.push_back(pDblE);
101         }
102         pDbl->killMe();
103         return types::Function::OK;
104     }
105
106     iMinRowsCols = std::min(pDbl->getRows(), pDbl->getCols());
107     pDblL = new types::Double(pDbl->getRows(), iMinRowsCols, pDbl->isComplex());
108     pDblU = new types::Double(iMinRowsCols, pDbl->getCols(), pDbl->isComplex());
109
110     if (pDbl->isComplex())
111     {
112         pdL = (double*) MALLOC(pDblL->getSize() * sizeof(doublecomplex));
113         pdU = (double*) MALLOC(pDblU->getSize() * sizeof(doublecomplex));
114     }
115     else
116     {
117         pdL = pDblL->get();
118         pdU = pDblU->get();
119     }
120
121     if (_iRetCount == 3)
122     {
123         pDblE = new types::Double(pDbl->getRows(), pDbl->getRows());
124     }
125
126     int iRet = iLuM(pData, pDbl->getRows(), pDbl->getCols(), pDbl->isComplex(), pdL, pdU, (pDblE ? pDblE->get() : NULL));
127
128     if (iRet != 0)
129     {
130         Scierror(999, _("%s: LAPACK error n°%d.\n"), "lu", iRet);
131         pDbl->killMe();
132         FREE((doublecomplex*)pdL);
133         FREE((doublecomplex*)pdU);
134         delete pDblL;
135         delete pDblU;
136         delete pDblE;
137         return types::Function::Error;
138     }
139
140     if (pDbl->isComplex())
141     {
142         vGetPointerFromDoubleComplex((doublecomplex*)pdL, pDblL->getSize(), pDblL->getReal(), pDblL->getImg());
143         FREE((doublecomplex*)pdL);
144         vGetPointerFromDoubleComplex((doublecomplex*)pdU, pDblU->getSize(), pDblU->getReal(), pDblU->getImg());
145         FREE((doublecomplex*)pdU);
146     }
147
148     if (pDbl->isComplex())
149     {
150         vFreeDoubleComplexFromPointer((doublecomplex*)pData);
151     }
152
153     pDbl->killMe();
154     out.push_back(pDblL);
155     out.push_back(pDblU);
156     if (_iRetCount == 3)
157     {
158         out.push_back(pDblE);
159     }
160
161     return types::Function::OK;
162 }
163 /*--------------------------------------------------------------------------*/
164