83eb39667a31dbdffd47490d995d0ee283f84083
[scilab.git] / scilab / modules / cacsd / sci_gateway / cpp / sci_arl2_ius.cpp
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2014 - Scilab Enterprises - 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 #include "cacsd_gw.hxx"
14 #include "function.hxx"
15 #include "double.hxx"
16 #include "polynom.hxx"
17 #include "string.hxx"
18
19 extern "C"
20 {
21 #include "sciprint.h"
22 #include "Scierror.h"
23 #include "localization.h"
24 #include "common_structure.h"
25 #include "elem_common.h"
26
27     extern void C2F(arl2)(double*, int*, double*, double*, int*, int*, double*, double*, int*, int*, int*, int*);
28     extern void C2F(arl2a)(double*, int*, double*, int*, int*, int*, int*, int*, int*, double*, int*);
29     extern void C2F(lq)(int*, double*, double*, double*, int*);
30     extern double C2F(phi)(double*, int*, double*, int*, double*);
31     extern void C2F(idegre)(double*, int*, int*);
32 }
33
34 /*--------------------------------------------------------------------------*/
35 types::Function::ReturnValue sci_arl2_ius(types::typed_list &in, int _iRetCount, types::typed_list &out)
36 {
37     types::Double* pDblY        = NULL;
38     types::Polynom* pPolyY      = NULL;
39     types::Polynom* pPolyDen    = NULL;
40     types::Double* pDblErr      = NULL;
41
42     double* pdblY   = NULL;
43     double* pdblDen = NULL;
44
45     int iOne        = 1;
46     int iVol1       = 0;
47     int iN          = 0;
48     bool bAll       = false;
49     int iRankDen    = 0;
50     double dErrl2   = 0;
51     int lunit       = 6;
52
53     C2F(arl2c).info = 0;
54
55     if (in.size() < 3 || in.size() > 5)
56     {
57         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "arl2_ius", 3, 5);
58         return types::Function::Error;
59     }
60
61     if (_iRetCount != 1 && _iRetCount > 3)
62     {
63         Scierror(78, _("%s: Wrong number of output argument(s): %d or %d expected.\n"), "arl2_ius", 1, 3);
64         return types::Function::Error;
65     }
66
67     /*** get inputs arguments ***/
68     // get Y
69     if (in[0]->isDouble())
70     {
71         pDblY = in[0]->getAs<types::Double>();
72         if (pDblY->isComplex())
73         {
74             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "arl2_ius", 1);
75             return types::Function::Error;
76         }
77
78         iVol1 = pDblY->getSize();
79         pdblY = pDblY->get();
80     }
81     else if (in[0]->isPoly())
82     {
83         pPolyY = in[0]->getAs<types::Polynom>();
84         if (pPolyY->isScalar() == false)
85         {
86             Scierror(999, _("%s: Wrong size for input argument #%d: A single polynom expected.\n"), "arl2_ius", 1);
87             return types::Function::Error;
88         }
89
90         if (pPolyY->isComplex())
91         {
92             Scierror(999, _("%s: Wrong type for input argument #%d: A real polynom expected.\n"), "arl2_ius", 1);
93             return types::Function::Error;
94         }
95
96         types::Double* pDblCoefY = pPolyY->get(0)->getCoef();
97         iVol1 = pDblCoefY->getSize();
98         pdblY = pDblCoefY->get();
99     }
100     else
101     {
102         Scierror(999, _("%s: Wrong type for input argument #%d: A Matrix or polynom expected.\n"), "arl2_ius", 1);
103         return types::Function::Error;
104     }
105
106     // get den <= useless in case "all" but it was does like that.
107     if (in[1]->isPoly() == false)
108     {
109         Scierror(999, _("%s: Wrong type for input argument #%d: A polynom expected.\n"), "arl2_ius", 2);
110         return types::Function::Error;
111     }
112
113     pPolyDen = in[1]->getAs<types::Polynom>();
114
115     if (pPolyDen->isScalar() == false)
116     {
117         Scierror(999, _("%s: Wrong size for input argument #%d: A single polynom expected.\n"), "arl2_ius", 2);
118         return types::Function::Error;
119     }
120
121     if (pPolyDen->isComplex())
122     {
123         Scierror(999, _("%s: Wrong value for input argument #%d: A real polynom expected.\n"), "arl2_ius", 2);
124         return types::Function::Error;
125     }
126
127     pPolyDen->getRealRank(&iRankDen);
128     pdblDen = pPolyDen->get(0)->getCoef()->get();
129     C2F(idegre)(pdblDen, &iRankDen, &iRankDen);
130     int iSize = iRankDen + 1;
131     double dblScal = 1.0 / pdblDen[iRankDen];
132     C2F(dscal)(&iSize, &dblScal, pdblDen, &iOne);
133
134     // get n
135     if (in[2]->isDouble() == false)
136     {
137         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "arl2_ius", 3);
138         return types::Function::Error;
139     }
140
141     types::Double* pDblN = in[2]->getAs<types::Double>();
142     iN = (int)pDblN->get(0);
143     if (iN < 1)
144     {
145         Scierror(999, _("%s: Wrong value for input argument #%d: More or equal to %d expected.\n"), "arl2_ius", 3, 1);
146         return types::Function::Error;
147     }
148
149     if (iN < iRankDen)
150     {
151         Scierror(999, _("%s: Wrong value for input argument #%d: More than degree of input argument #%d expected.\n"), "arl2_ius", 3, 2);
152         return types::Function::Error;
153     }
154
155     // get "all" and/or imp
156     if (in.size() == 4)
157     {
158         if (in[3]->isString()) // get "all"
159         {
160             types::String* pStrAll = in[3]->getAs<types::String>();
161             if (pStrAll->isScalar() == false)
162             {
163                 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar string expected.\n"), "arl2_ius", 4);
164                 return types::Function::Error;
165             }
166
167             if (wcscmp(pStrAll->get(0), L"all") != 0)
168             {
169                 Scierror(999, _("%s: Wrong value for input argument #%d: 'all' expected.\n"), "arl2_ius", 4);
170                 return types::Function::Error;
171             }
172
173             bAll = true;
174         }
175         else if (in[3]->isDouble()) // get imp
176         {
177             types::Double* pDblImp = in[3]->getAs<types::Double>();
178             C2F(arl2c).info = (int)pDblImp->get(0);
179             if (C2F(arl2c).info < 0)
180             {
181                 Scierror(999, _("%s: Wrong value for input argument #%d: Positive value expected.\n"), "arl2_ius", 4);
182                 return types::Function::Error;
183             }
184         }
185         else
186         {
187             Scierror(999, _("%s: Wrong type for input argument #%d: A scalar or string expected.\n"), "arl2_ius", 4);
188             return types::Function::Error;
189         }
190     }
191     else if (in.size() == 5)
192     {
193         // get imp
194         if (in[3]->isDouble() == false)
195         {
196             Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "arl2_ius", 4);
197             return types::Function::Error;
198         }
199
200         types::Double* pDblImp = in[3]->getAs<types::Double>();
201         C2F(arl2c).info = (int)pDblImp->get(0);
202         if (C2F(arl2c).info < 0)
203         {
204             Scierror(999, _("%s: Wrong value for input argument #%d: Positive value expected.\n"), "arl2_ius", 4);
205             return types::Function::Error;
206         }
207
208         // get "all"
209         if (in[4]->isString() == false)
210         {
211             Scierror(999, _("%s: Wrong type for input argument #%d: A string expected.\n"), "arl2_ius", 5);
212             return types::Function::Error;
213         }
214
215         types::String* pStrAll = in[4]->getAs<types::String>();
216         if (pStrAll->isScalar() == false)
217         {
218             Scierror(999, _("%s: Wrong size for input argument #%d: A scalar string expected.\n"), "arl2_ius", 5);
219             return types::Function::Error;
220         }
221
222         if (wcscmp(pStrAll->get(0), L"all") != 0)
223         {
224             Scierror(999, _("%s: Wrong value for input argument #%d: 'all' expected.\n"), "arl2_ius", 5);
225             return types::Function::Error;
226         }
227
228         bAll = true;
229     }
230
231     /*** perform operations ***/
232     int iNg = iVol1 - 1;
233     if (bAll)
234     {
235         // look for "all" solutions
236         int iNsol = 0;
237         int iMxsol = 20;
238         double* pdblDenTemp = new double[iMxsol * (iN + 1)];
239         int iWorkSize = 34 + 34 * iN + 7 * iNg + iN * iNg + iN * iN * (iNg + 2) + 4 * (iN + 1) * iMxsol;
240         double* pdblWork = new double[iWorkSize];
241         iWorkSize = 29 + iN * iN + 4 * iN + 2 * iMxsol;
242         int* piWork = new int[iWorkSize];
243
244         C2F(arl2a)(pdblY, &iVol1, pdblDenTemp, &iMxsol, &iNsol, &iN,
245                    &(C2F(arl2c).info), &(C2F(arl2c).ierr), &lunit, pdblWork, piWork);
246
247         delete[] pdblWork;
248         delete[] piWork;
249
250         if (C2F(arl2c).ierr)
251         {
252             if (C2F(arl2c).ierr == 3)
253             {
254                 Scierror(999, _("%s: Loop on two orders detected.\n"), "arl2a");
255                 return types::Function::Error;
256             }
257             else if (C2F(arl2c).ierr == 4)
258             {
259                 Scierror(999, _("%s: Impossible to reach required order.\n"), "arl2a");
260                 return types::Function::Error;
261             }
262             else if (C2F(arl2c).ierr == 5)
263             {
264                 Scierror(999, _("%s: Failure when looking for the intersection with domains bounds.\n"), "arl2a");
265                 return types::Function::Error;
266             }
267             else if (C2F(arl2c).ierr == 7)
268             {
269                 Scierror(999, _("%s: Too many solutions found.\n"), "arl2a");
270                 return types::Function::Error;
271             }
272         }
273
274         /*** retrun output arguments ***/
275         // retrun denominators
276         double** pdblAllCoeff = new double*[iNsol];
277         int iRank = iN + 1;
278         int* piRank = new int[iNsol];
279         for (int i = 0; i < iNsol; i++)
280         {
281             piRank[i] = iRank;
282         }
283
284         types::Polynom* pPolyDenOut = new types::Polynom(pPolyDen->getVariableName(), iNsol, 1, piRank);
285         for (int i = 0; i < iNsol; i++)
286         {
287             double* pdblDenOut = pPolyDenOut->get(i)->getCoef()->get();
288             C2F(dcopy)(&iN, pdblDenTemp + i, &iMxsol, pdblDenOut, &iOne);
289             pdblDenOut[iN] = 1;
290             pdblAllCoeff[i] = pdblDenOut;
291         }
292
293         delete[] pdblDenTemp;
294
295         out.push_back(pPolyDenOut);
296
297         // retrun numerators
298         if (_iRetCount > 1)
299         {
300             for (int i = 0; i < iNsol; i++)
301             {
302                 piRank[i] = iN;
303             }
304
305             C2F(no2f).gnrm = sqrt(C2F(no2f).gnrm);
306
307             types::Polynom* pPolyNumOut = new types::Polynom(pPolyDen->getVariableName(), iNsol, 1, piRank);
308             for (int i = 0; i < iNsol; i++)
309             {
310                 double* pdblNumOut = pPolyNumOut->get(i)->getCoef()->get();
311                 double* pdblWork = new double[iN + iNg + 1];
312                 C2F(lq)(&iN, pdblAllCoeff[i], pdblWork, pdblY, &iNg);
313                 C2F(dscal)(&iN, &(C2F(no2f).gnrm), pdblWork, &iOne);
314                 C2F(dcopy)(&iN, pdblWork, &iOne, pdblNumOut, &iOne);
315                 delete[] pdblWork;
316             }
317
318             out.push_back(pPolyNumOut);
319         }
320
321         // return error
322         if (_iRetCount > 2)
323         {
324             pDblErr = new types::Double(iNsol, 1);
325             double* pdblerr = pDblErr->get();
326             double* pdblWork = new double[iN + iNg + 1];
327             for (int i = 0; i < iNsol; i++)
328             {
329                 pdblerr[i] = sqrt(C2F(phi)(pdblAllCoeff[i], &iN, pdblY, &iNg, pdblWork)) * C2F(no2f).gnrm;
330             }
331
332             delete[] pdblWork;
333
334             out.push_back(pDblErr);
335         }
336
337         delete[] piRank;
338         delete[] pdblAllCoeff;
339     }
340     else
341     {
342         // look for a solution
343         int iSizeNum = (std::max)(iN, iRankDen);
344         double* pdblNum = new double[iSizeNum];
345         int iWorkSize = 32 + 32 * iN + 7 * iNg + iN * iNg + iN * iN * (iNg + 2);
346         double* pdblWork = new double[iWorkSize];
347         iWorkSize = 29 + iN * iN + 4 * iN;
348         int* piWork = new int[iWorkSize];
349
350         int iSizeTemp = (std::max)(iRankDen, iN) + 1;
351         double* pDblDenTemp = new double[iSizeTemp];
352         memset(pDblDenTemp, 0x00, iSizeTemp * sizeof(double));
353         C2F(dcopy)(&iSize, pdblDen, &iOne, pDblDenTemp, &iOne);
354
355         C2F(arl2)(pdblY, &iVol1, pdblNum, pDblDenTemp, &iRankDen, &iN, &dErrl2, pdblWork, piWork,
356                   &(C2F(arl2c).info), &(C2F(arl2c).ierr), &lunit);
357
358
359         delete[] pdblWork;
360         delete[] piWork;
361
362         if (C2F(arl2c).ierr != 0)
363         {
364             if (C2F(arl2c).ierr == 3)
365             {
366                 sciprint(_("%s: Loop on two orders detected.\n"), "arl2");
367             }
368             else if (C2F(arl2c).ierr == 4)
369             {
370                 sciprint(_("%s: Impossible to reach required order.\n   previous order computed solution returned.\n"), "arl2");
371             }
372             else if (C2F(arl2c).ierr == 5)
373             {
374                 sciprint(_("%s: Failure when looking for the intersection with domains boundaries.\n   previous order computed solution returned.\n"), "arl2");
375             }
376             else if (C2F(arl2c).ierr == 7)
377             {
378                 Scierror(999, _("%s: too many solutions found\n"), "arl2");
379                 delete[] pdblNum;
380                 delete[] pDblDenTemp;
381                 return types::Function::Error;
382             }
383         }
384
385         /*** retrun output arguments ***/
386         // retrun denominator
387         int iRank = iN + 1;
388         types::Polynom* pPolyDenOut = new types::Polynom(pPolyDen->getVariableName(), 1, 1, &iRank);
389         double* pdblDenOut = pPolyDenOut->get(0)->getCoef()->get();
390         C2F(dcopy)(&iRank, pDblDenTemp, &iOne, pdblDenOut, &iOne);
391         out.push_back(pPolyDenOut);
392
393         // retrun numerator
394         if (_iRetCount > 1)
395         {
396             types::Polynom* pPolyNumOut = new types::Polynom(pPolyDen->getVariableName(), 1, 1, &iN);
397             double* pdblNumOut = pPolyNumOut->get(0)->getCoef()->get();
398             C2F(dcopy)(&iN, pdblNum, &iOne, pdblNumOut, &iOne);
399             out.push_back(pPolyNumOut);
400         }
401
402         // return error
403         if (_iRetCount > 2)
404         {
405             out.push_back(new types::Double(dErrl2));
406         }
407
408         delete[] pDblDenTemp;
409         delete[] pdblNum;
410     }
411
412     return types::Function::OK;
413 }
414 /*--------------------------------------------------------------------------*/