fd75b416ba274d043ea0e2279b18de9395003bfd
[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                 int iSize1 = piRank[0][i] + 1;
149                 int iSize2 = piRank[1][i] + 1;
150                 int iSize3 = piRank[2][i] + 1;
151                 C2F(residu)(pdblInR[0][i], &iSize1, pdblInR[1][i], &iSize2,
152                             pdblInR[2][i], &iSize3, &v, &dblEps, &iErr);
153                 if (iErr)
154                 {
155                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "residu");
156                     throw iErr;
157                 }
158
159                 pdblOut[i] = v;
160             }
161         }
162         else
163         {
164             // complex case
165             pDblOut = new types::Double(iRows[0], iCols[0], true);
166             double* pdblOutR = pDblOut->get();
167             double* pdblOutI = pDblOut->getImg();
168
169             for (int i = 0; i < 3; i++)
170             {
171                 if (pdblInI[i] == NULL)
172                 {
173                     pdblInI[i] = new double*[iSize];
174                     for (int j = 0; j < iSize; j++)
175                     {
176                         int iLen = piRank[i][j] + 1;
177                         pdblInI[i][j] = new double[iLen];
178                         memset(pdblInI[i][j], 0x00, iLen * sizeof(double));
179                     }
180
181                     isDeletable[i] = true;
182                 }
183             }
184
185             for (int i = 0; i < iRows[0] * iCols[0]; i++)
186             {
187                 int iErr    = 0;
188                 double real = 0;
189                 double imag = 0;
190
191                 C2F(wesidu)(pdblInR[0][i], pdblInI[0][i], (piRank[0]) + i,
192                             pdblInR[1][i], pdblInI[1][i], (piRank[1]) + i,
193                             pdblInR[2][i], pdblInI[2][i], (piRank[2]) + i,
194                             &real, &imag, &dblEps, &iErr);
195
196                 if (iErr)
197                 {
198                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "wesidu");
199                     throw iErr;
200                 }
201
202                 pdblOutR[i] = real;
203                 pdblOutI[i] = imag;
204             }
205         }
206     }
207     catch (int error)
208     {
209         iError = error;
210     }
211
212     // free memory
213     for (int i = 0; i < 3; i++)
214     {
215         if (pDblIn[i])
216         {
217             delete pDblIn[i];
218         }
219
220         if (pPoly[i])
221         {
222             delete pPoly[i];
223         }
224
225         if (piRank[i])
226         {
227             delete[] piRank[i];
228         }
229
230         if (pdblInR[i])
231         {
232             delete[] pdblInR[i];
233         }
234
235         if (isDeletable[i])
236         {
237             for (int j = 0; j < iSize; j++)
238             {
239                 delete[] pdblInI[i][j];
240             }
241         }
242
243         if (pdblInI[i])
244         {
245             delete[] pdblInI[i];
246         }
247     }
248
249     /*** retrun output arguments ***/
250     if (iError)
251     {
252         return types::Function::Error;
253     }
254
255     out.push_back(pDblOut);
256     return types::Function::OK;
257 }
258 /*--------------------------------------------------------------------------*/