linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / src / cpp / schurSelect.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
4 *
5 * This file must be used under the terms of the CeCILL.
6 * This source file is licensed as described in the file COPYING, which
7 * you should have received as part of this distribution.  The terms
8 * are also available at
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12 /*--------------------------------------------------------------------------*/
13
14 #include "configvariable.hxx"
15 #include "callable.hxx"
16 #include "runvisitor.hxx"
17
18 extern "C"
19 {
20     #include "scischur.h"
21     #include "schurSelect.h"
22 }
23
24 /*--------------------------------------------------------------------------*/
25 int schurSelect(types::Double** _pDblIn, types::Double** _pDblOut, bool _bIsComplex, bool _bIsDiscrete, bool _bIsContinu, ConfigVariable::EntryPointStr* pStrFunction)
26 {
27     int info                    = 0;
28     int* pBwork                 = NULL;
29     int iWorksize               = 0;
30     double* pRwork              = NULL;
31     doublecomplex* pCplxWork    = NULL;
32     int iDim                    = 0;
33     int iCols                   = _pDblIn[0]->getCols();
34     types::Callable* pCall      = ConfigVariable::getSchurFunction();
35
36     doublecomplex* pDataInDoublecomplex[2]     = {NULL, NULL};
37     doublecomplex* pDataOutDoublecomplex[2]    = {NULL, NULL};
38
39     pBwork = (int*)MALLOC((_pDblIn[1] ? 2 : 1) * iCols * sizeof(int));
40     if(pBwork == NULL)
41     {
42         return -1;
43     }
44
45     const char* jobL = _pDblOut[0]  ? "V":"N";
46     const char* jobR = _pDblOut[1]  ? "V":"N";
47     const char* sort = (pStrFunction || pCall || _bIsDiscrete || _bIsContinu) ? "S":"N";
48
49     if(_pDblIn[1] == NULL && _bIsComplex == false)
50     {//dgees
51         double* pWR = (double*)MALLOC(iCols * sizeof(double));
52         double* pWI = (double*)MALLOC(iCols * sizeof(double));
53         pRwork = allocDgeesWorkspace(iCols, &iWorksize);
54
55         if(pWR == NULL || pWI == NULL || pRwork == NULL)
56         {
57             return -1;
58         }
59
60         if(_bIsDiscrete)
61         {
62             C2F(dgees)(jobL, sort, schur_sb02mw, &iCols, _pDblIn[0]->getReal(), &iCols, &iDim, pWR, pWI, _pDblOut[0]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
63         }
64         else if(_bIsContinu)
65         {
66             C2F(dgees)(jobL, sort, schur_sb02mv, &iCols, _pDblIn[0]->getReal(), &iCols, &iDim, pWR, pWI, _pDblOut[0]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
67         }
68         else if(pCall)
69         {
70             C2F(dgees)(jobL, sort, schur_dgees, &iCols, _pDblIn[0]->getReal(), &iCols, &iDim, pWR, pWI, _pDblOut[0]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
71         }
72         else if(pStrFunction)
73         {
74             C2F(dgees)(jobL, sort, (schur_dgees_t)pStrFunction->functionPtr, &iCols, _pDblIn[0]->getReal(), &iCols, &iDim, pWR, pWI, _pDblOut[0]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
75         }
76         else
77         {
78             C2F(dgees)(jobL, sort, NULL, &iCols, _pDblIn[0]->getReal(), &iCols, &iDim, pWR, pWI, _pDblOut[0]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
79         }
80
81         if(_pDblOut[2])
82         {
83             _pDblOut[2]->set(0, (double)iDim);
84         }
85
86             FREE(pWR);
87             FREE(pWI);
88         FREE(pRwork);
89     }
90
91     if(_pDblIn[1] == NULL && _bIsComplex)
92     {//zgees
93         doublecomplex* pW = NULL;
94         pRwork      = (double*)MALLOC(iCols * sizeof(double));
95         pW          = (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex));
96         pCplxWork   = allocZgeesWorkspace(iCols, &iWorksize);
97
98         if(pRwork == NULL || pW == NULL || pCplxWork == NULL)
99         {
100             return -1;
101         }
102
103         pDataInDoublecomplex[0] = oGetDoubleComplexFromPointer(_pDblIn[0]->getReal(), _pDblIn[0]->getImg(), _pDblIn[0]->getSize());
104         pDataOutDoublecomplex[0] = oGetDoubleComplexFromPointer(_pDblOut[0]->getReal(), _pDblOut[0]->getImg(), _pDblOut[0]->getSize());
105
106         if(_bIsDiscrete)
107         {
108             C2F(zgees)(jobL, sort, schur_zb02mw, &iCols, pDataInDoublecomplex[0], &iCols, &iDim, pW, pDataOutDoublecomplex[0], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
109         }
110         else if(_bIsContinu)
111         {
112             C2F(zgees)(jobL, sort, schur_zb02mv, &iCols, pDataInDoublecomplex[0], &iCols, &iDim, pW, pDataOutDoublecomplex[0], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
113         }
114         else if(pCall)
115         {
116             C2F(zgees)(jobL, sort, schur_zgees, &iCols, pDataInDoublecomplex[0], &iCols, &iDim, pW, pDataOutDoublecomplex[0], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
117         }
118         else if(pStrFunction)
119         {
120             C2F(zgees)(jobL, sort, (schur_zgees_t)pStrFunction->functionPtr, &iCols, pDataInDoublecomplex[0], &iCols, &iDim, pW, pDataOutDoublecomplex[0], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
121         }
122         else
123         {
124             C2F(zgees)(jobL, sort, NULL, &iCols, pDataInDoublecomplex[0], &iCols, &iDim, pW, pDataOutDoublecomplex[0], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
125         }
126
127         if(_pDblOut[2])
128         {
129             _pDblOut[2]->set(0, (double)iDim);
130         }
131
132         vGetPointerFromDoubleComplex(pDataInDoublecomplex[0], _pDblIn[0]->getSize(), _pDblIn[0]->getReal(), _pDblIn[0]->getImg());
133         vGetPointerFromDoubleComplex(pDataOutDoublecomplex[0], _pDblOut[0]->getSize(), _pDblOut[0]->getReal(), _pDblOut[0]->getImg());
134
135             FREE(pW);
136         FREE(pRwork);
137         FREE(pCplxWork);
138         vFreeDoubleComplexFromPointer(pDataInDoublecomplex[0]);
139         vFreeDoubleComplexFromPointer(pDataOutDoublecomplex[0]);
140     }
141
142     if(_pDblIn[1] && _bIsComplex == false)
143     {//dgges
144             double* pAlphaR = (double*)MALLOC(iCols * sizeof(double));
145             double* pAlphaI = (double*)MALLOC(iCols * sizeof(double));
146             double* pBeta   = (double*)MALLOC(iCols * sizeof(double));
147         pRwork = allocDggesWorkspace(iCols, &iWorksize);
148
149         if(pAlphaR == NULL || pAlphaI == NULL || pBeta == NULL || pRwork == NULL)
150         {
151             return -1;
152         }
153
154         if(_bIsDiscrete)
155         {
156             C2F(dgges)(jobL, jobR, sort, schur_sb02ox, &iCols, _pDblIn[0]->getReal(), &iCols, _pDblIn[1]->getReal(), &iCols, &iDim, pAlphaR, pAlphaI, pBeta, _pDblOut[0]->get(), &iCols, _pDblOut[1]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
157         }
158         else if(_bIsContinu)
159         {
160             C2F(dgges)(jobL, jobR, sort, schur_sb02ow, &iCols, _pDblIn[0]->getReal(), &iCols, _pDblIn[1]->getReal(), &iCols, &iDim, pAlphaR, pAlphaI, pBeta, _pDblOut[0]->get(), &iCols, _pDblOut[1]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
161         }
162         else if(pCall)
163         {
164             C2F(dgges)(jobL, jobR, sort, schur_dgges, &iCols, _pDblIn[0]->getReal(), &iCols, _pDblIn[1]->getReal(), &iCols, &iDim, pAlphaR, pAlphaI, pBeta, _pDblOut[0]->get(), &iCols, _pDblOut[1]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
165         }
166         else if(pStrFunction)
167         {
168             C2F(dgges)(jobL, jobR, sort, (schur_dgges_t)pStrFunction->functionPtr, &iCols, _pDblIn[0]->getReal(), &iCols, _pDblIn[1]->getReal(), &iCols, &iDim, pAlphaR, pAlphaI, pBeta, _pDblOut[0]->get(), &iCols, _pDblOut[1]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
169         }
170         else
171         {
172             C2F(dgges)(jobL, jobR, sort, NULL, &iCols, _pDblIn[0]->getReal(), &iCols, _pDblIn[1]->getReal(), &iCols, &iDim, pAlphaR, pAlphaI, pBeta, _pDblOut[0]->get(), &iCols, _pDblOut[1]->get(), &iCols, pRwork, &iWorksize, pBwork, &info);
173         }
174
175         if(_pDblOut[2])
176         {
177             _pDblOut[2]->set(0, (double)iDim);
178         }
179
180             FREE(pAlphaR);
181             FREE(pAlphaI);
182             FREE(pBeta);
183         FREE(pRwork);
184     }
185
186     if(_pDblIn[1] && _bIsComplex)
187     {//zgges
188             doublecomplex* pAlpha   = (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex));
189             doublecomplex* pBeta    = (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex));
190         pRwork                  = (double*) MALLOC(8 * iCols * sizeof(double)); 
191         pCplxWork               = allocZggesWorkspace(iCols, &iWorksize);
192
193         if(pRwork == NULL || pAlpha == NULL || pBeta == NULL || pCplxWork == NULL)
194         {
195             return -1;
196         }
197
198         pDataInDoublecomplex[0] = oGetDoubleComplexFromPointer(_pDblIn[0]->getReal(), _pDblIn[0]->getImg(), _pDblIn[0]->getSize());
199         pDataOutDoublecomplex[0] = oGetDoubleComplexFromPointer(_pDblOut[0]->getReal(), _pDblOut[0]->getImg(), _pDblOut[0]->getSize());
200         pDataInDoublecomplex[1] = oGetDoubleComplexFromPointer(_pDblIn[1]->getReal(), _pDblIn[1]->getImg(), _pDblIn[1]->getSize());
201         pDataOutDoublecomplex[1] = oGetDoubleComplexFromPointer(_pDblOut[1]->getReal(), _pDblOut[1]->getImg(), _pDblOut[1]->getSize());
202
203         if(_bIsDiscrete)
204         {
205             C2F(zgges)(jobL, jobR, sort, schur_zb02ox, &iCols, pDataInDoublecomplex[0], &iCols, pDataInDoublecomplex[1], &iCols, &iDim, pAlpha, pBeta, pDataOutDoublecomplex[0], &iCols, pDataOutDoublecomplex[1], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
206         }
207         else if(_bIsContinu)
208         {
209             C2F(zgges)(jobL, jobR, sort, schur_zb02ow, &iCols, pDataInDoublecomplex[0], &iCols, pDataInDoublecomplex[1], &iCols, &iDim, pAlpha, pBeta, pDataOutDoublecomplex[0], &iCols, pDataOutDoublecomplex[1], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
210         }
211         else if(pCall)
212         {
213             C2F(zgges)(jobL, jobR, sort, schur_zgges, &iCols, pDataInDoublecomplex[0], &iCols, pDataInDoublecomplex[1], &iCols, &iDim, pAlpha, pBeta, pDataOutDoublecomplex[0], &iCols, pDataOutDoublecomplex[1], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
214         }
215         else if(pStrFunction)
216         {
217             C2F(zgges)(jobL, jobR, sort, (schur_zgges_t)pStrFunction->functionPtr, &iCols, pDataInDoublecomplex[0], &iCols, pDataInDoublecomplex[1], &iCols, &iDim, pAlpha, pBeta, pDataOutDoublecomplex[0], &iCols, pDataOutDoublecomplex[1], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
218         }
219         else
220         {
221             C2F(zgges)(jobL, jobR, sort, NULL, &iCols, pDataInDoublecomplex[0], &iCols, pDataInDoublecomplex[1], &iCols, &iDim, pAlpha, pBeta, pDataOutDoublecomplex[0], &iCols, pDataOutDoublecomplex[1], &iCols, pCplxWork, &iWorksize, pRwork, pBwork, &info);
222         }
223
224         if(_pDblOut[2])
225         {
226             _pDblOut[2]->set(0, (double)iDim);
227         }
228
229         vGetPointerFromDoubleComplex(pDataInDoublecomplex[0], _pDblIn[0]->getSize(), _pDblIn[0]->getReal(), _pDblIn[0]->getImg());
230         vGetPointerFromDoubleComplex(pDataOutDoublecomplex[0], _pDblOut[0]->getSize(), _pDblOut[0]->getReal(), _pDblOut[0]->getImg());
231         vGetPointerFromDoubleComplex(pDataInDoublecomplex[1], _pDblIn[1]->getSize(), _pDblIn[1]->getReal(), _pDblIn[1]->getImg());
232         vGetPointerFromDoubleComplex(pDataOutDoublecomplex[1], _pDblOut[1]->getSize(), _pDblOut[1]->getReal(), _pDblOut[1]->getImg());
233
234         FREE(pRwork);
235         FREE(pCplxWork);
236         vFreeDoubleComplexFromPointer(pDataInDoublecomplex[0]);
237         vFreeDoubleComplexFromPointer(pDataOutDoublecomplex[0]);
238         vFreeDoubleComplexFromPointer(pDataInDoublecomplex[1]);
239         vFreeDoubleComplexFromPointer(pDataOutDoublecomplex[1]);
240     }
241
242     FREE(pBwork);
243     return info;
244 }
245 /*--------------------------------------------------------------------------*/
246 double* allocDgeesWorkspace(int iCols, int* allocated)
247 {
248     int info;
249     int query = -1;
250     double optim;
251     double* ret = NULL;
252
253     C2F(dgees)("V", "N", NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, NULL, &iCols, &optim, &query, NULL, &info);
254
255     *allocated = (int)optim;
256     ret = (double*) MALLOC(*allocated * sizeof(double));
257
258     if(!ret)
259     {
260         *allocated = 3 * iCols;
261         ret = (double*) MALLOC(*allocated * sizeof(double));
262
263         if(!ret)
264         {
265             *allocated = 0;
266         }
267     }
268     return ret;
269 }
270
271 doublecomplex* allocZgeesWorkspace(int iCols, int* allocated)
272 {
273     int info;
274     int query = -1;
275     doublecomplex optim;
276     doublecomplex* ret = NULL;
277
278     C2F(zgees)("V", "N", NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, &optim, &query, NULL, NULL, &info);
279
280     *allocated = (int)optim.r;
281     ret = (doublecomplex*) MALLOC(*allocated * sizeof(doublecomplex));
282
283     if(!ret)
284     {
285         *allocated = 2 * iCols;
286         ret = (doublecomplex*) MALLOC(*allocated * sizeof(doublecomplex));
287         if(!ret)
288         {
289             *allocated = 0;
290         }
291     }
292     return ret;
293 }
294
295 double* allocDggesWorkspace(int iCols, int* allocated)
296 {
297     int info;
298     int query = -1;
299     double optim;
300     double* ret = NULL;
301
302     C2F(dgges)("V", "V", "N", NULL, &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, NULL, NULL, &iCols, NULL, &iCols, &optim, &query, NULL, &info);
303
304     *allocated = (int)optim;
305     ret = (double*)MALLOC(*allocated * sizeof(double));
306
307     if(!ret)
308     {
309         *allocated = 8 * iCols + 16;
310         ret = (double*)MALLOC(*allocated * sizeof(double));
311
312         if(!ret)
313         {
314             *allocated = 0;
315         }
316     }
317     return ret;
318 }
319
320 doublecomplex* allocZggesWorkspace(int iCols, int* allocated)
321 {
322     int info;
323     int query = -1;
324     doublecomplex optim;
325     doublecomplex* ret = NULL;
326
327     C2F(zgges)("V", "V", "N", NULL, &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, NULL, &iCols, NULL, &iCols, &optim, &query, NULL, NULL, &info);
328
329     *allocated = (int)optim.r;
330     ret = (doublecomplex*) MALLOC(*allocated * sizeof(doublecomplex));
331
332     if(!ret)
333     {
334         *allocated = 2 * iCols;
335         ret = (doublecomplex*) MALLOC(*allocated * sizeof(doublecomplex));
336         if(!ret)
337         {
338             *allocated = 0;
339         }
340     }
341     return ret;
342 }
343 /*--------------------------------------------------------------------------*/