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