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