linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_schur.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 #include "configvariable.hxx"
21 #include "callable.hxx"
22
23 extern "C"
24 {
25 #include "localization.h"
26 #include "Scierror.h"
27 #include "scischur.h"
28 #include "schurSelect.h"
29 }
30
31 /*--------------------------------------------------------------------------*/
32 types::Function::ReturnValue sci_schur(types::typed_list &in, int _iRetCount, types::typed_list &out)
33 {
34     types::Double* pDbl[2]          = {NULL, NULL};
35     types::String* pStr             = NULL;
36     types::Double* pDblOut[3]       = {NULL, NULL, NULL};
37     types::Callable* pFunction      = NULL;
38     ConfigVariable::EntryPointStr* pStrFunction = NULL;
39
40     int  iCase                      = 1;
41     bool bComplexArgs               = false;
42     bool bIsRealStr                 = false;
43     bool bIsComplexStr              = false;
44     bool bIsContinuStr              = false;
45     bool bIsDiscreteStr             = false;
46
47 // *** check number of input args. ***
48     if(in.size() < 1 && in.size() > 3)
49     {
50         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d to %d expected.\n"), L"schur", 1, 3);
51         return types::Function::Error;
52     }
53
54 // *** check type of input args and get it. ***
55     if(in[0]->isDouble() == false)
56     {
57         std::wstring wstFuncName = L"%"  + in[0]->getShortTypeStr() + L"_schur";
58         return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
59     }
60
61     pDbl[0] = in[0]->getAs<types::Double>();
62
63     if(pDbl[0]->getRows() != pDbl[0]->getCols()) // square matrix
64     {
65         ScierrorW(20, _W("%ls: Wrong type for argument %d: Square matrix expected.\n"), L"schur", 1);
66         return types::Function::Error;
67     }
68
69     if((pDbl[0]->getRows() == -1) || (pDbl[0]->getCols() == -1)) // manage eye case
70     {
71         ScierrorW(271,_W("%ls: Size varying argument a*eye(), (arg %d) not allowed here.\n"), L"schur", 1);
72         return types::Function::Error;
73     }
74
75     pDbl[0] = pDbl[0]->clone()->getAs<types::Double>(); // pDbl will be modified
76
77     if(in.size() > 1)
78     {
79         if(in[1]->isDouble())
80         {
81             pDbl[1] = in[1]->getAs<types::Double>();
82
83             if(pDbl[1]->getRows() != pDbl[1]->getCols()) // square matrix
84             {
85                 ScierrorW(20, _W("%ls: Wrong type for argument %d: Square matrix expected.\n"), L"schur", 2);
86                 return types::Function::Error;
87             }
88
89             if((pDbl[1]->getRows() == -1) || (pDbl[1]->getCols() == -1)) // manage eye case
90             {
91                 ScierrorW(271,_W("%ls: Size varying argument a*eye(), (arg %d) not allowed here.\n"), L"schur", 2);
92                 return types::Function::Error;
93             }
94
95             if(pDbl[0]->getCols() != pDbl[1]->getCols())
96             {
97                 ScierrorW(267,_W("%ls: Arg %d and arg %d must have equal dimensions.\n"), L"schur", 1, 2);
98                 return types::Function::Error;
99             }
100
101             pDbl[1] = pDbl[1]->clone()->getAs<types::Double>(); // pDbl will be modified
102
103             iCase = 11;
104         }
105         else if(in[1]->isString())
106         {
107             pStr = in[1]->getAs<types::String>();
108             bIsRealStr = !wcscmp(pStr->get(0), L"r") || !wcscmp(pStr->get(0), L"real");
109             bIsComplexStr = !wcscmp(pStr->get(0), L"comp") || !wcscmp(pStr->get(0), L"complex");
110
111             if(bIsComplexStr)
112             {
113                 pDbl[0]->setComplex(true); // pDbl[0] is a clone of in[0]
114             }
115
116             if(bIsRealStr || bIsComplexStr)
117             {
118                 iCase = 1;
119             }
120             else
121             {
122                 iCase = 12;
123             }
124         }
125         else if(in[1]->isCallable())
126         {
127             pFunction = in[1]->getAs<types::Callable>();
128             iCase = 12;
129         }
130         else
131         {
132             std::wstring wstFuncName = L"%"  + in[1]->getShortTypeStr() + L"_schur";
133             return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
134         }
135     }
136
137     if(in.size() == 3)
138     {
139         if(in[2]->isString() == false && in[2]->isCallable() == false)
140         {
141             std::wstring wstFuncName = L"%"  + in[2]->getShortTypeStr() + L"_schur";
142             return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
143         }
144
145         if(in[2]->isString())
146         {
147             pStr = in[2]->getAs<types::String>();
148         }
149         else
150         {
151             pFunction = in[2]->getAs<types::Callable>();
152         }
153
154         iCase = 112;
155     }
156
157 // *** check number of output args according the input args. ***
158
159     // iCase represents the type of args input.
160     // 1 = double, 2 = string, so 112 = double double string.
161     switch(iCase)
162     {
163         case 1: // double
164         {
165             if(_iRetCount > 2)
166             {
167                 ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d to %d expected.\n"), L"schur", 1, 2);
168                 return types::Function::Error;
169             }
170             break;
171         }
172         case 11: // double double
173         {
174             if(_iRetCount != 2 && _iRetCount != 4)
175             {
176                 ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d or %d expected.\n"), L"schur", 2, 4);
177                 return types::Function::Error;
178             }
179             break;
180         }
181         case 12: // double string
182         {
183             if(_iRetCount < 2 && _iRetCount > 3)
184             {
185                 ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d to %d expected.\n"), L"schur", 2, 3);
186                 return types::Function::Error;
187             }
188             break;
189         }
190         default://case 112: // double double string
191         {
192             if(_iRetCount > 4) // in doc, 5 output args are possible ?!?
193             {
194                 ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d to %d expected.\n"), L"schur", 1, 4);
195                 return types::Function::Error;
196             }
197             break;
198         }
199     }
200
201 // *** check empty matrix ***
202
203     // iCase represents the type of args input.
204     // 1 = double, 2 = string, so 112 = double double string.
205     switch(iCase)
206     {
207         case 11: // double double
208         case 1: // double
209         {
210             if(pDbl[0]->getCols() == 0)
211             {
212                 for(int i=0; i<_iRetCount; i++)
213                 {
214                     out.push_back(types::Double::Empty());
215                 }
216                 return types::Function::OK;
217             }
218             break;
219         }
220         case 12: // double string
221         {
222             if(pDbl[0]->getCols() == 0)
223             {
224                 types::Double* zero = new types::Double(0);
225
226                 for(int i = 0; i < _iRetCount; i++)
227                 {
228                     if(i == 1 && !bIsComplexStr && !bIsRealStr)
229                     {
230                         out.push_back(zero);
231                     }
232                     else
233                     {
234                         out.push_back(types::Double::Empty());
235                     }
236                 }
237                 return types::Function::OK;
238             }
239             break;
240         }
241         case 112: // double double string
242         {
243             if(pDbl[0]->getCols() == 0)
244             {
245                 types::Double* zero = new types::Double(0);
246
247                 for(int i=1; i<_iRetCount; i++)
248                 {
249                     out.push_back(types::Double::Empty());
250                 }
251
252                 if(_iRetCount > 1)
253                 {
254                     out.push_back(zero);
255                 }
256                 else
257                 {
258                     out.push_back(types::Double::Empty());
259                 }
260
261                 return types::Function::OK;
262             }
263             break;
264         }
265     }
266
267 // *** perform operations. ***
268
269     bComplexArgs = pDbl[0]->isComplex() || (pDbl[1] && pDbl[1]->isComplex()) || bIsComplexStr;
270
271     if(bIsRealStr && bComplexArgs)
272     {
273         ScierrorW(202,_W("%ls: Wrong type for input argument #%d: Real matrix expected.\n"), L"schur", 1);
274         return types::Function::Error;
275     }
276
277     for(int i = 0; i < 2; i++)
278     {
279         pDblOut[i] = new types::Double(pDbl[0]->getCols(), pDbl[0]->getCols(), bComplexArgs);
280     }
281
282     // iCase represents the type of args input.
283     // 1 = double, 2 = string, so 112 = double double string.
284     switch(iCase)
285     {
286         case 112: // double double string || double double function
287         case 12: // double string || double function
288         {
289             if(pStr)
290             {
291                 bIsContinuStr = !wcscmp(pStr->get(0), L"c") ||
292                                 !wcscmp(pStr->get(0), L"cont") ||
293                                 !wcscmp(pStr->get(0), L"zb02ow") || // two matrix, complex case
294                                 !wcscmp(pStr->get(0), L"zb02mv") || // one matrix, complex case
295                                 !wcscmp(pStr->get(0), L"sb02ow") || // two matrix, real case
296                                 !wcscmp(pStr->get(0), L"sb02mv");   // one matrix, real case
297
298                 bIsDiscreteStr = !wcscmp(pStr->get(0), L"d") ||
299                                 !wcscmp(pStr->get(0), L"disc") ||
300                                 !wcscmp(pStr->get(0), L"zb02ox") || // two matrix, complex case
301                                 !wcscmp(pStr->get(0), L"zb02mw") || // one matrix, complex case
302                                 !wcscmp(pStr->get(0), L"sb02ox") || // two matrix, real case
303                                 !wcscmp(pStr->get(0), L"sb02mw");   // one matrix, real case
304
305                 if(bIsContinuStr == false && bIsDiscreteStr == false)
306                 {
307                     pStrFunction = ConfigVariable::getEntryPoint(pStr->get(0));
308                     if(pStrFunction == NULL)
309                     {
310                         ScierrorW(999,_W("%ls: Subroutine not found: %ls\n"), L"schur", pStr->get(0));
311                         return types::Function::Error;
312                     }
313                 }
314             }
315             else if(pFunction)
316             {
317                 ConfigVariable::setSchurFunction(pFunction);
318             }
319
320             pDblOut[2] = new types::Double(1,1); // Dim
321             break;
322         }
323         default:// case 1 and 11
324             break;
325     }
326
327     int iRet = schurSelect(pDbl, pDblOut, bComplexArgs, bIsDiscreteStr, bIsContinuStr, pStrFunction);
328
329     if(iRet)
330     {
331         printf("failed !!!  %d",iRet);
332         return types::Function::Error;
333     }
334
335     // return result(s)
336     switch(iCase)
337     {
338         case 1: // double
339         {
340             if(_iRetCount == 2)
341             {
342                 out.push_back(pDblOut[0]);
343             }
344
345             out.push_back(pDbl[0]); // pDbl[0] has been overwritten by its real Schur form T.
346             break;
347         }
348         case 11: // double double
349         {
350             for(int i = 0; i < 2; i++)
351             {
352                 out.push_back(pDbl[i]);
353             }
354
355             if(_iRetCount == 4)
356             {
357                 out.push_back(pDblOut[0]);
358                 if(_iRetCount > 1)
359                 {
360                     out.push_back(pDblOut[1]);
361                 }
362             }
363
364             break;
365         }
366         case 12: // double string || double function
367         {
368             if(_iRetCount < 2)
369             {
370                 out.push_back(pDbl[0]);
371             }
372             else
373             {
374                 out.push_back(pDblOut[0]);
375                 out.push_back(pDblOut[2]);
376
377                 if(_iRetCount == 3)
378                 {
379                     out.push_back(pDbl[0]);
380                 }
381             }
382             break;
383         }
384         case 112: // double double string || double double function
385         {
386             switch(_iRetCount)
387             {
388                 case 0 :
389                 case 1 : break; // dim
390                 case 4 : // As Es Z dim
391                 {
392                     for(int i = 0; i < 2; i++)
393                     {
394                         out.push_back(pDbl[i]);
395                     }
396                     out.push_back(pDblOut[1]);
397                     break;
398                 }
399                 case 3 : // Q Z dim
400                 {
401                     out.push_back(pDblOut[0]);
402                 }
403                 case 2 : // Z dim
404                 {
405                     out.push_back(pDblOut[1]);
406                     break;
407                 }
408             }
409             out.push_back(pDblOut[2]);
410             break;
411         }
412     }
413
414     ConfigVariable::setSchurFunction(NULL);
415     return types::Function::OK;
416 }
417 /*--------------------------------------------------------------------------*/