refactoring polynom.
[scilab.git] / scilab / modules / cacsd / sci_gateway / cpp / sci_rtitr.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
18 extern "C"
19 {
20 #include "sciprint.h"
21 #include "Scierror.h"
22 #include "localization.h"
23 #include "dmp2pm.h"
24
25     extern void C2F(rtitr)( int*, int*, int*, double*, int*, int*, double*, int*, int*, double*,
26                             double*, int*, double*, double*, int*, int*, int*, double*, int*);
27
28 }
29
30 /*--------------------------------------------------------------------------*/
31 types::Function::ReturnValue sci_rtitr(types::typed_list &in, int _iRetCount, types::typed_list &out)
32 {
33     types::Double* pDblUp = NULL;
34     types::Double* pDblYp = NULL;
35
36     double* pdblUp      = NULL;
37     double* pdblYp      = NULL;
38     double* pdblU       = NULL;
39     double** pdblNum    = NULL;
40     double** pdblDen    = NULL;
41     int* piRankDen      = NULL;
42     int* piRankNum      = NULL;
43     int iMaxRankNum     = 0;
44     int iMaxRankDen     = 0;
45
46     int iRowsDen = 0;
47     int iSizeDen = 0;
48     int iRowsNum = 0;
49     int iColsNum = 0;
50     int iSizeNum = 0;
51     int iRowsU   = 0;
52     int iColsU   = 0;
53     int iOne     = 1;
54     int iJob     = 1;
55     int iErr     = 0;
56
57     if (in.size() != 3 && in.size() != 5)
58     {
59         Scierror(77, _("%s: Wrong number of input argument(s): %d or %d expected.\n"), "rtitr", 3, 5);
60         return types::Function::Error;
61     }
62
63     if (_iRetCount < 1)
64     {
65         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "rtitr", 1);
66         return types::Function::Error;
67     }
68
69     /*** get inputs arguments ***/
70     // get up and yp
71     if (in.size() == 5)
72     {
73         if (in[4]->isDouble() == false)
74         {
75             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "rtitr", 5);
76             return types::Function::Error;
77         }
78
79         pDblYp = in[4]->getAs<types::Double>();
80
81         if (pDblYp->isComplex())
82         {
83             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "rtitr", 5);
84             return types::Function::Error;
85         }
86
87         pdblYp = pDblYp->get();
88
89         if (in[3]->isDouble() == false)
90         {
91             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "rtitr", 4);
92             return types::Function::Error;
93         }
94
95         pDblUp = in[3]->getAs<types::Double>();
96
97         if (pDblUp->isComplex())
98         {
99             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "rtitr", 4);
100             return types::Function::Error;
101         }
102
103         pdblUp = pDblUp->get();
104
105         iJob = 2;
106     }
107
108     // get Num
109     if (in[0]->isDouble())
110     {
111         types::Double* pDblNum = in[0]->getAs<types::Double>();
112         if (pDblNum->isComplex())
113         {
114             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "rtitr", 1);
115             return types::Function::Error;
116         }
117
118         iSizeNum    = pDblNum->getSize();
119         iRowsNum    = pDblNum->getRows();
120         iColsNum    = pDblNum->getCols();
121         iMaxRankNum = 0;
122         pdblNum     = new double*[iSizeNum];
123         double* pdbl = pDblNum->get();
124
125         for (int i = 0; i < iSizeNum; i++)
126         {
127             pdblNum[i] = pdbl + i;
128         }
129     }
130     else if (in[0]->isPoly())
131     {
132         types::Polynom* pPolyNum = in[0]->getAs<types::Polynom>();
133         if (pPolyNum->isComplex())
134         {
135             Scierror(999, _("%s: Wrong type for input argument #%d: A real polynom expected.\n"), "rtitr", 1);
136             return types::Function::Error;
137         }
138
139         iSizeNum    = pPolyNum->getSize();
140         iRowsNum    = pPolyNum->getRows();
141         iColsNum    = pPolyNum->getCols();
142         iMaxRankNum = pPolyNum->getMaxRank();
143         piRankNum   = new int[iSizeNum];
144         pPolyNum->getRank(piRankNum);
145
146         pdblNum     = new double*[iSizeNum];
147         for (int i = 0; i < iSizeNum; i++)
148         {
149             pdblNum[i] = pPolyNum->get(i)->get();
150         }
151     }
152     else
153     {
154         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or polynom expected.\n"), "rtitr", 1);
155         return types::Function::Error;
156     }
157
158     // get Den
159     if (in[1]->isDouble())
160     {
161         types::Double* pDblDen = in[1]->getAs<types::Double>();
162         if (pDblDen->isComplex())
163         {
164             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "rtitr", 2);
165             return types::Function::Error;
166         }
167
168         if (pDblDen->getRows() != pDblDen->getCols())
169         {
170             Scierror(999, _("%s: Wrong size for input argument #%d: A square matrix expected.\n"), "rtitr", 2);
171             return types::Function::Error;
172         }
173
174         iSizeDen    = pDblDen->getSize();
175         iRowsDen    = pDblDen->getRows();
176         iMaxRankDen = 0;
177         pdblDen     = new double*[iSizeDen];
178         double* pdbl = pDblDen->get();
179
180         for (int i = 0; i < iSizeDen; i++)
181         {
182             pdblDen[i] = pdbl + i;
183         }
184     }
185     else if (in[1]->isPoly())
186     {
187         types::Polynom* pPolyDen = in[1]->getAs<types::Polynom>();
188         if (pPolyDen->isComplex())
189         {
190             Scierror(999, _("%s: Wrong type for input argument #%d: A real polynom expected.\n"), "rtitr", 2);
191             return types::Function::Error;
192         }
193
194         if (pPolyDen->getRows() != pPolyDen->getCols())
195         {
196             Scierror(999, _("%s: Wrong size for input argument #%d: A square matrix expected.\n"), "rtitr", 2);
197             return types::Function::Error;
198         }
199
200         iSizeDen    = pPolyDen->getSize();
201         iRowsDen    = pPolyDen->getRows();
202         iMaxRankDen = pPolyDen->getMaxRank();
203         piRankDen   = new int[iSizeDen];
204         pPolyDen->getRank(piRankDen);
205
206         pdblDen     = new double*[iSizeDen];
207         for (int i = 0; i < iSizeDen; i++)
208         {
209             pdblDen[i] = pPolyDen->get(i)->get();
210         }
211     }
212     else
213     {
214         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or polynom expected.\n"), "rtitr", 2);
215         return types::Function::Error;
216     }
217
218     // get u
219     if (in[2]->isDouble() == false)
220     {
221         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "rtitr", 3);
222         return types::Function::Error;
223     }
224
225     types::Double* pDblU = in[2]->getAs<types::Double>();
226
227     if (pDblU->isComplex())
228     {
229         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "rtitr", 3);
230         return types::Function::Error;
231     }
232
233     iRowsU = pDblU->getRows();
234     iColsU = pDblU->getCols();
235
236     if (iRowsDen != iRowsNum || iColsNum != iRowsU)
237     {
238         Scierror(60, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "rtitr");
239         return types::Function::Error;
240     }
241
242     pdblU = pDblU->get();
243
244     // check size of input argument 4 and 5
245     if (in.size() == 5)
246     {
247         if (pDblYp->getRows() != iRowsDen && pDblYp->getRows() != 0)
248         {
249             Scierror(60, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "rtitr");
250             return types::Function::Error;
251         }
252
253         if (pDblYp->getCols() != iMaxRankDen)
254         {
255             Scierror(60, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "rtitr");
256             return types::Function::Error;
257         }
258
259         if (pDblUp->getRows() != iColsNum && pDblUp->getRows() != 0)
260         {
261             Scierror(60, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "rtitr");
262             return types::Function::Error;
263         }
264
265         if (pDblUp->getCols() != iMaxRankDen)
266         {
267             Scierror(60, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "rtitr");
268             return types::Function::Error;
269         }
270     }
271
272     /*** perform operations ***/
273     // output
274     int iRowsOut            = iColsU + iMaxRankDen - iMaxRankNum;
275     types::Double* pDblOut  = new types::Double(iRowsDen, iRowsOut);
276     double* pdblOut         = pDblOut->get();
277
278     // working array
279     double* pdblWork    = new double[iRowsDen];
280     int* piWork         = new int[iRowsDen];
281
282     // converted a matrix of polynom to a polynomial matrix
283     double* pdblDenMat = dmp2pm(pdblDen, iSizeDen, piRankDen, iMaxRankDen);
284     double* pdblNumMat = dmp2pm(pdblNum, iSizeNum, piRankNum, iMaxRankNum);
285
286     // perform operations
287     C2F(rtitr)(&iRowsU, &iRowsDen, &iColsU, pdblNumMat, &iRowsDen, &iMaxRankNum, pdblDenMat, &iRowsDen, &iMaxRankDen,
288                pdblUp, pdblU, &iRowsU, pdblYp, pdblOut, &iRowsDen, &iJob, piWork, pdblWork, &iErr);
289
290     // check error
291     if (iErr)
292     {
293         if (iErr == 1)
294         {
295             char strValue[256];
296             sprintf(strValue, "%lf", pdblWork[0]);
297             sciprint(_("Warning :\n"));
298             sciprint(_("matrix is close to singular or badly scaled. rcond = %s\n"), strValue);
299             iErr = 0;
300         }
301         else if (iErr == 2)
302         {
303             Scierror(19, _("Problem is singular.\n"), "rtitr");
304         }
305     }
306
307     // free memory
308     free(pdblDenMat);
309     free(pdblNumMat);
310     delete[] pdblWork;
311     delete[] piWork;
312     delete[] pdblDen;
313     delete[] pdblNum;
314
315     if (piRankDen)
316     {
317         delete[] piRankDen;
318     }
319
320     if (piRankNum)
321     {
322         delete[] piRankNum;
323     }
324
325     /*** retrun output arguments ***/
326     if (iErr)
327     {
328         return types::Function::Error;
329     }
330
331     out.push_back(pDblOut);
332     return types::Function::OK;
333 }
334 /*--------------------------------------------------------------------------*/