refactoring polynom.
[scilab.git] / scilab / modules / cacsd / sci_gateway / cpp / sci_residu.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 "Scierror.h"
21 #include "localization.h"
22 #include "elem_common.h"
23
24     extern void C2F(residu)(double*, int*, double*, int*, double*, int*, double*, double*, int*);
25     extern void C2F(wesidu)(double*, double*, int*, double*, double*, int*,
26                             double*, double*, int*, double*, double*, double*, int*);
27 }
28
29 /*--------------------------------------------------------------------------*/
30 types::Function::ReturnValue sci_residu(types::typed_list &in, int _iRetCount, types::typed_list &out)
31 {
32     int iRows[3]             = {0, 0, 0};
33     int iCols[3]             = {0, 0, 0};
34     int iComplex[3]          = {0, 0, 0};
35     int* piRank[3]           = {NULL, NULL, NULL};
36     double** pdblInR[3]      = {NULL, NULL, NULL};
37     double** pdblInI[3]      = {NULL, NULL, NULL};
38     bool isDeletable[3]      = {false, false, false};
39     types::Double* pDblIn[3] = {NULL, NULL, NULL};
40     types::Polynom* pPoly[3] = {NULL, NULL, NULL};
41     types::Double* pDblOut   = NULL;
42
43     char cP         = 'p';
44     double dblEps   = C2F(dlamch)(&cP, 1L);
45     double dZero    = 0;
46     int iOne        = 1;
47     int iSize       = 0;
48     int iError      = 0;
49
50     if (in.size() != 3)
51     {
52         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "residu", 3);
53         return types::Function::Error;
54     }
55
56     if (_iRetCount > 1)
57     {
58         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "residu", 1);
59         return types::Function::Error;
60     }
61
62     try
63     {
64         /*** get inputs arguments ***/
65         for (int i = 0; i < 3; i++)
66         {
67             if (in[i]->isDouble())
68             {
69                 pDblIn[i] = in[i]->clone()->getAs<types::Double>();
70                 iRows[i] = pDblIn[i]->getRows();
71                 iCols[i] = pDblIn[i]->getCols();
72
73                 iSize = pDblIn[i]->getSize();
74                 piRank[i] = new int[iSize];
75                 memset(piRank[i], 0x00, iSize * sizeof(int));
76
77                 pdblInR[i]  = new double*[iSize];
78                 double* pdbl = pDblIn[i]->get();
79                 for (int j = 0; j < iSize; i++)
80                 {
81                     pdblInR[i][j] = pdbl + j;
82                 }
83
84                 if (pDblIn[i]->isComplex())
85                 {
86                     pdblInI[i]  = new double*[iSize];
87                     double* pdbl = pDblIn[i]->get();
88                     for (int j = 0; j < iSize; i++)
89                     {
90                         pdblInI[i][j] = pdbl + j;
91                     }
92                 }
93             }
94             else if (in[i]->isPoly())
95             {
96                 pPoly[i] = in[i]->clone()->getAs<types::Polynom>();
97                 iRows[i] = pPoly[i]->getRows();
98                 iCols[i] = pPoly[i]->getCols();
99
100                 iSize = pPoly[i]->getSize();
101                 piRank[i] = new int[iSize];
102                 pPoly[i]->getRank(piRank[i]);
103
104                 pdblInR[i] = new double*[iSize];
105                 if (pPoly[i]->isComplex())
106                 {
107                     pdblInI[i] = new double*[iSize];
108                     for (int j = 0; j < iSize; j++)
109                     {
110                         pdblInR[i][j] = pPoly[i]->get(j)->get();
111                         pdblInI[i][j] = pPoly[i]->get(j)->getImg();
112                     }
113                 }
114                 else
115                 {
116                     for (int j = 0; j < iSize; j++)
117                     {
118                         pdblInR[i][j] = pPoly[i]->get(j)->get();
119                     }
120                 }
121             }
122             else
123             {
124                 Scierror(999, _("%s: Wrong type for input argument #%d: A Matrix or polynom expected.\n"), "residu", i + 1);
125                 throw 1;
126             }
127         }
128
129         if (iRows[0] != iRows[1] || iCols[0] != iCols[1] || iRows[0] != iRows[2] || iCols[0] != iCols[2])
130         {
131             Scierror(999, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "residu");
132             throw 1;
133         }
134
135         /*** perform operations ***/
136         if (pdblInI[0] == NULL && pdblInI[1] == NULL && pdblInI[2] == NULL)
137         {
138             // real case
139             pDblOut = new types::Double(iRows[0], iCols[0]);
140             double* pdblOut = pDblOut->get();
141             for (int i = 0; i < iRows[0] * iCols[0]; i++)
142             {
143                 int iErr = 0;
144                 double v = 0;
145
146                 int iSize1 = piRank[0][i] + 1;
147                 int iSize2 = piRank[1][i] + 1;
148                 int iSize3 = piRank[2][i] + 1;
149                 C2F(residu)(pdblInR[0][i], &iSize1, pdblInR[1][i], &iSize2,
150                             pdblInR[2][i], &iSize3, &v, &dblEps, &iErr);
151                 if (iErr)
152                 {
153                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "residu");
154                     throw iErr;
155                 }
156
157                 pdblOut[i] = v;
158             }
159         }
160         else
161         {
162             // complex case
163             pDblOut = new types::Double(iRows[0], iCols[0], true);
164             double* pdblOutR = pDblOut->get();
165             double* pdblOutI = pDblOut->getImg();
166
167             for (int i = 0; i < 3; i++)
168             {
169                 if (pdblInI[i] == NULL)
170                 {
171                     pdblInI[i] = new double*[iSize];
172                     for (int j = 0; j < iSize; j++)
173                     {
174                         int iLen = piRank[i][j] + 1;
175                         pdblInI[i][j] = new double[iLen];
176                         memset(pdblInI[i][j], 0x00, iLen * sizeof(double));
177                     }
178
179                     isDeletable[i] = true;
180                 }
181             }
182
183             for (int i = 0; i < iRows[0] * iCols[0]; i++)
184             {
185                 int iErr    = 0;
186                 double real = 0;
187                 double imag = 0;
188
189                 C2F(wesidu)(pdblInR[0][i], pdblInI[0][i], (piRank[0]) + i,
190                             pdblInR[1][i], pdblInI[1][i], (piRank[1]) + i,
191                             pdblInR[2][i], pdblInI[2][i], (piRank[2]) + i,
192                             &real, &imag, &dblEps, &iErr);
193
194                 if (iErr)
195                 {
196                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "wesidu");
197                     throw iErr;
198                 }
199
200                 pdblOutR[i] = real;
201                 pdblOutI[i] = imag;
202             }
203         }
204     }
205     catch (int error)
206     {
207         iError = error;
208     }
209
210     // free memory
211     for (int i = 0; i < 3; i++)
212     {
213         if (pDblIn[i])
214         {
215             delete pDblIn[i];
216         }
217
218         if (pPoly[i])
219         {
220             delete pPoly[i];
221         }
222
223         if (piRank[i])
224         {
225             delete[] piRank[i];
226         }
227
228         if (pdblInR[i])
229         {
230             delete[] pdblInR[i];
231         }
232
233         if (isDeletable[i])
234         {
235             for (int j = 0; j < iSize; j++)
236             {
237                 delete[] pdblInI[i][j];
238             }
239         }
240
241         if (pdblInI[i])
242         {
243             delete[] pdblInI[i];
244         }
245     }
246
247     /*** retrun output arguments ***/
248     if (iError)
249     {
250         return types::Function::Error;
251     }
252
253     out.push_back(pDblOut);
254     return types::Function::OK;
255 }
256 /*--------------------------------------------------------------------------*/