memory leak fixed in svd gateway.
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_svd.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 "svd.h"
26 #include "doublecomplex.h"
27 }
28 /*--------------------------------------------------------------------------*/
29 /*
30   s=svd(X): [R: min(rows, cols) x 1]
31   [U,S,V]=svd(X) [U,S,V]:  [ [C|R: rows x rows], [R: rows x cols ],  [R|C: cols x cols ] ]
32   [U,S,V]=svd(X,0) (obsolete) [U,S,V]=svd(X,"e"): [ [C|R:rows x min(rows,cols)], [R: min(rows,cols) x min(rows,cols)], [C|R:cols x min(rows,cols)] ]
33   [U,S,V,rk]=svd(X [,tol]) : cf. supra, rk[R 1 x 1]
34
35   /!\ Contrary to specifications (@ http://www.scilab.org/product/man/svd.html )
36   , previous version was accepting Lhs==2. Worse : tests were using this undocumented behavior.
37   implementation and tests have been fixed according to the specification.
38
39 */
40 extern "C"
41 {
42     extern int C2F(vfinite)(int *n, double *v);
43 }
44 /*--------------------------------------------------------------------------*/
45
46 types::Function::ReturnValue sci_svd(types::typed_list &in, int _iRetCount, types::typed_list &out)
47 {
48     types::Double* pDbl     = NULL;
49     types::Double* ptrsU    = NULL;
50     types::Double* ptrsV    = NULL;
51     types::Double* ptrS     = NULL;
52     types::Double* pRk      = NULL;
53     types::Double* pSV      = NULL;
54     double* pData           = NULL;
55     int economy             = 0;
56     int totalsize           = 0;
57     int iRet                = 0;
58     double tol              = 0.;
59
60     if (in.size() != 1 && in.size() != 2)
61     {
62         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "svd", 1, 2);
63         return types::Function::Error;
64     }
65
66     if (_iRetCount > 4)
67     {
68         Scierror(78, _("%s: Wrong number of output argument(s): At least %d expected.\n"), "svd", 4);
69         return types::Function::Error;
70     }
71
72     if (in[0]->isDouble() == false)
73     {
74         ast::ExecVisitor exec;
75         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_svd";
76         return Overload::call(wstFuncName, in, _iRetCount, out, &exec);
77     }
78     pDbl = in[0]->clone()->getAs<types::Double>();
79
80     if (in.size() == 2)
81     {
82         if (in[1]->isString() == true)
83         {
84             if (_iRetCount == 4)
85             {
86                 Scierror(78, _("%s: Wrong number of output argument(s): %d or %d expected.\n"), "svd", 1, 3);
87                 return types::Function::Error;
88             }
89             types::String* pStr = in[1]->getAs<types::String>();
90             if ((wcslen(pStr->get(0)) == 1) && (pStr->get(0)[0] == L'e'))
91             {
92                 economy = 1;
93             }
94         }
95         if (in[1]->isDouble() == true)
96         {
97             /* no further testing for "old Economy size:  [U,S,V]=svd(A,0) " */
98             if (_iRetCount == 3)
99             {
100                 economy = 1;
101             }
102         }
103     }
104
105     if (pDbl->isEmpty()) /* empty matrix */
106     {
107         for (int i = 0; i < _iRetCount - 1; i++)
108         {
109             out.push_back(types::Double::Empty());
110         }
111
112         if (_iRetCount == 4)
113         {
114             types::Double* pDZero = new types::Double(1, 1);
115             pDZero->set(0, 0);
116             out.push_back(pDZero);
117         }
118         else
119         {
120             out.push_back(types::Double::Empty());
121         }
122
123         delete pDbl;
124         return types::Function::OK;
125     }
126
127     if ((pDbl->getRows() == -1) || (pDbl->getCols() == -1)) // manage eye case
128     {
129         Scierror(271, _("%s: Size varying argument a*eye(), (arg %d) not allowed here.\n"), "svd", 1);
130         delete pDbl;
131         return types::Function::Error;
132     }
133
134     if (pDbl->isComplex())
135     {
136         pData = (double *)oGetDoubleComplexFromPointer(pDbl->getReal(), pDbl->getImg(), pDbl->getSize());
137     }
138     else
139     {
140         pData = pDbl->get();
141     }
142
143     totalsize = pDbl->getSize() * (pDbl->isComplex() ? 2 : 1);
144     if (C2F(vfinite)(&totalsize, pData) == false)
145     {
146         Scierror(264, _("%s: Wrong value for argument %d: Must not contain NaN or Inf.\n"), "svd", 1);
147         delete pDbl;
148         return types::Function::Error;
149     }
150
151     switch (_iRetCount)
152     {
153         case 0:
154         case 1:
155         {
156             pSV = new types::Double(std::min(pDbl->getRows(), pDbl->getCols()), 1, false);
157             iRet = iSvdM(pData, pDbl->getRows(), pDbl->getCols(), pDbl->isComplex(), economy, tol, pSV->get(), NULL, NULL, NULL, NULL);
158         }
159         break;
160         case 4:
161         {
162             if ((in.size() == 2) && (in[1]->isDouble()))
163             {
164                 tol = in[1]->getAs<types::Double>()->get(0);
165             }
166             pRk = new types::Double(1, 1, false);
167         }
168         case 2:
169         case 3:
170         {
171             int economyRows = economy ? std::min(pDbl->getRows(), pDbl->getCols()) : pDbl->getRows();
172             int economyCols = economy ? std::min(pDbl->getRows(), pDbl->getCols()) : pDbl->getCols();
173
174             ptrsU = new types::Double(pDbl->getRows(), economyRows, pDbl->isComplex());
175             ptrS  = new types::Double(economyRows, economyCols, false);
176             ptrsV = new types::Double(pDbl->getCols(), economyCols, pDbl->isComplex());
177
178             if (pDbl->isComplex())
179             {
180                 double* U = (double*)MALLOC(pDbl->getRows() * economyRows * sizeof(doublecomplex));
181                 double* V = (double*)MALLOC(pDbl->getCols() * economyCols * sizeof(doublecomplex));
182
183                 iRet = iSvdM(pData, pDbl->getRows(), pDbl->getCols(), true /*isComplex*/, economy, tol, NULL, U, ptrS->get(), V, (pRk ? pRk->get() : NULL));
184
185                 vGetPointerFromDoubleComplex((doublecomplex*)U, ptrsU->getSize(), ptrsU->getReal(), ptrsU->getImg());
186                 vFreeDoubleComplexFromPointer((doublecomplex*)U);
187                 vGetPointerFromDoubleComplex((doublecomplex*)V, ptrsV->getSize(), ptrsV->getReal(), ptrsV->getImg());
188                 vFreeDoubleComplexFromPointer((doublecomplex*)V);
189             }
190             else
191             {
192                 iRet = iSvdM(pData, pDbl->getRows(), pDbl->getCols(), false /*isComplex*/, economy, tol, NULL, ptrsU->get(), ptrS->get(), ptrsV->get(), (pRk ? pRk->get() : NULL));
193             }
194         }
195         break;
196         // default: // makes at the beginning of this gateway
197     }
198
199     if (iRet != 0)
200     {
201         if (iRet == -1)
202         {
203             Scierror(999, _("%s: Cannot allocate more memory.\n"), "svd");
204         }
205         else
206         {
207             Scierror(24, _("%s: Convergence problem...\n"), "svd");
208         }
209
210         delete pDbl;
211         return types::Function::Error;
212     }
213
214     if (pDbl->isComplex())
215     {
216         vFreeDoubleComplexFromPointer((doublecomplex*)pData);
217     }
218
219     switch (_iRetCount)
220     {
221         case 4:
222         {
223             out.push_back(ptrsU);
224             out.push_back(ptrS);
225             out.push_back(ptrsV);
226             out.push_back(pRk);
227             break;
228         }
229         case 3:
230         {
231             out.push_back(ptrsU);
232             out.push_back(ptrS);
233             out.push_back(ptrsV);
234             break;
235         }
236         case 2:
237         {
238             out.push_back(ptrsU);
239             out.push_back(ptrS);
240             delete ptrsV;
241             break;
242         }
243         case 1:
244             out.push_back(pSV);
245             // default: // makes at the beginning of this gateway
246     }
247
248     delete pDbl;
249     return types::Function::OK;
250 }
251 /*--------------------------------------------------------------------------*/
252