* Bug 15248: lsq() was leaking memory
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_lsq.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 "lsq.h"
27 #include "doublecomplex.h"
28 }
29 /*--------------------------------------------------------------------------*/
30
31 types::Function::ReturnValue sci_lsq(types::typed_list &in, int _iRetCount, types::typed_list &out)
32 {
33     types::Double* pDbl[2]      = {NULL, NULL};
34     types::Double* pDblResult   = NULL;
35     double* pData[2]            = {NULL, NULL};
36     double* pResult             = NULL;
37     double* pdTol               = NULL;
38     bool bComplexArgs           = false;
39     int iRank                   = 0;
40     double dblTol               = 0.0;
41
42     if (in.size() < 2 || in.size() > 3)
43     {
44         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "lsq", 2, 3);
45         return types::Function::Error;
46     }
47
48     if (_iRetCount > 2)
49     {
50         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "lsq", 1, 2);
51         return types::Function::Error;
52     }
53
54     if ((in[0]->isDouble() == false))
55     {
56         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lsq";
57         return Overload::call(wstFuncName, in, _iRetCount, out);
58     }
59
60     if (in.size() <=  3)
61     {
62         if ((in[1]->isDouble() == false))
63         {
64             std::wstring wstFuncName = L"%" + in[1]->getShortTypeStr() + L"_lsq";
65             return Overload::call(wstFuncName, in, _iRetCount, out);
66         }
67     }
68
69     if (in.size() == 3)
70     {
71         if ((in[2]->isDouble() == false) || (in[2]->getAs<types::Double>()->isComplex()) || (in[2]->getAs<types::Double>()->isScalar() == false))
72         {
73             Scierror(256, _("%s: Wrong type for input argument #%d: A Real expected.\n"), "lsq", 3);
74             return types::Function::Error;
75         }
76
77         dblTol = in[2]->getAs<types::Double>()->get(0);
78         pdTol = &dblTol;
79     }
80
81     pDbl[0] = in[0]->getAs<types::Double>();
82     pDbl[1] = in[1]->getAs<types::Double>();
83
84     if (pDbl[0]->getRows() != pDbl[1]->getRows())
85     {
86         Scierror(265, _("%s: %s and %s must have equal number of rows.\n"), "lsq", "A", "B");
87         return types::Function::Error;
88     }
89
90     if ((pDbl[0]->getCols() == 0) || (pDbl[1]->getCols() == 0))
91     {
92         out.push_back(types::Double::Empty());
93         if (_iRetCount == 2)
94         {
95             out.push_back(types::Double::Empty());
96         }
97         return types::Function::OK;
98     }
99
100     if (pDbl[0]->isComplex() || pDbl[1]->isComplex())
101     {
102         bComplexArgs = true;
103     }
104
105     pDbl[0] = pDbl[0]->clone();
106     pDbl[1] = pDbl[1]->clone();
107
108     for (int i = 0; i < 2; i++)
109     {
110         if (pDbl[i]->getCols() == -1)
111         {
112             for (int j=0; j<=i; j++)
113             {
114                 pDbl[j]->killMe();
115             }
116             Scierror(271, _("%s: Size varying argument a*eye(), (arg %d) not allowed here.\n"), "lsq", i + 1);
117             return types::Function::Error;
118         }
119
120         if (bComplexArgs)
121         {
122             pData[i] = (double*)oGetDoubleComplexFromPointer(pDbl[i]->getReal(), pDbl[i]->getImg(), pDbl[i]->getSize());
123             if (!pData[i])
124             {
125                 for (int j=0; j<=i; j++)
126                 {
127                     pDbl[j]->killMe();
128                 }
129                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "lsq");
130                 return types::Function::Error;
131             }
132         }
133         else
134         {
135             pData[i] = pDbl[i]->getReal();
136         }
137     }
138
139     pDblResult = new types::Double(pDbl[0]->getCols(), pDbl[1]->getCols(), bComplexArgs);
140
141     if (bComplexArgs)
142     {
143         pResult = (double*)MALLOC(pDbl[0]->getCols() * pDbl[1]->getCols() * sizeof(doublecomplex));
144     }
145     else
146     {
147         pResult = pDblResult->get();
148     }
149
150     int iRet = iLsqM(pData[0], pDbl[0]->getRows(), pDbl[0]->getCols(), pData[1], pDbl[1]->getCols(), bComplexArgs, pResult, pdTol, ((_iRetCount == 2) ? &iRank : NULL));
151
152     pDbl[0]->killMe();
153     pDbl[1]->killMe();
154
155     if (iRet != 0)
156     {
157         if (iRet == -1)
158         {
159             Scierror(999, _("%s: Allocation failed.\n"),  "lsq");
160         }
161         else
162         {
163             Scierror(999, _("%s: LAPACK error n°%d.\n"),  "lsq", iRet);
164         }
165
166         if (bComplexArgs)
167         {
168             vFreeDoubleComplexFromPointer((doublecomplex*)pResult);
169             vFreeDoubleComplexFromPointer((doublecomplex*)pData[0]);
170             vFreeDoubleComplexFromPointer((doublecomplex*)pData[1]);
171         }
172         pDblResult->killMe();
173         return types::Function::Error;
174     }
175
176     if (bComplexArgs)
177     {
178         vGetPointerFromDoubleComplex((doublecomplex*)(pResult), pDblResult->getSize(), pDblResult->getReal(), pDblResult->getImg());
179         vFreeDoubleComplexFromPointer((doublecomplex*)pResult);
180         vFreeDoubleComplexFromPointer((doublecomplex*)pData[0]);
181         vFreeDoubleComplexFromPointer((doublecomplex*)pData[1]);
182     }
183
184     out.push_back(pDblResult);
185     if (_iRetCount == 2)
186     {
187         types::Double* pDblRank = new types::Double(1, 1);
188         pDblRank->set(0, iRank);
189         out.push_back(pDblRank);
190     }
191
192     return types::Function::OK;
193 }
194 /*--------------------------------------------------------------------------*/