fc1f25e8b16d909d54958e7ab7b43d404d5ffc4e
[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]->getRealRank(piRank[i]);
103
104                 pdblInR[i] = new double*[iSize];
105                 for(int j = 0; j < iSize; j++)
106                 {
107                     pdblInR[i][j] = pPoly[i]->get(j)->getCoef()->get();
108                 }
109
110                 if (pPoly[i]->isComplex())
111                 {
112                     pdblInI[i] = new double*[iSize];
113                     for(int j = 0; j < iSize; j++)
114                     {
115                         pdblInI[i][j] = pPoly[i]->get(j)->getCoef()->getImg();
116                     }
117                 }
118             }
119             else
120             {
121                 Scierror(999, _("%s: Wrong type for input argument #%d: A Matrix or polynom expected.\n"), "residu", i+1);
122                 throw 1;
123             }
124         }
125
126         if(iRows[0] != iRows[1] || iCols[0] != iCols[1] || iRows[0] != iRows[2] || iCols[0] != iCols[2])
127         {
128             Scierror(999, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "residu");
129             throw 1;
130         }
131
132         /*** perform operations ***/
133         if(pdblInI[0] == NULL && pdblInI[1] == NULL && pdblInI[2] == NULL)
134         {// real case
135             pDblOut = new types::Double(iRows[0], iCols[0]);
136             double* pdblOut = pDblOut->get();
137             for(int i = 0; i < iRows[0] * iCols[0]; i++)
138             {
139                 int iErr = 0;
140                 double v = 0;
141
142                 int iSize1 = piRank[0][i] + 1;
143                 int iSize2 = piRank[1][i] + 1;
144                 int iSize3 = piRank[2][i] + 1;
145                 C2F(residu)(pdblInR[0][i], &iSize1, pdblInR[1][i], &iSize2,
146                             pdblInR[2][i], &iSize3, &v, &dblEps, &iErr);
147                 if(iErr)
148                 {
149                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "residu");
150                     throw iErr;
151                 }
152
153                 pdblOut[i] = v;
154             }
155         }
156         else
157         {// complex case
158             pDblOut = new types::Double(iRows[0], iCols[0], true);
159             double* pdblOutR = pDblOut->get();
160             double* pdblOutI = pDblOut->getImg();
161
162             for(int i= 0; i < 3; i++)
163             {
164                 if(pdblInI[i] == NULL)
165                 {
166                     pdblInI[i] = new double*[iSize];
167                     for(int j = 0; j < iSize; j++)
168                     {
169                         int iLen = piRank[i][j] + 1;
170                         pdblInI[i][j] = new double[iLen];
171                         memset(pdblInI[i][j], 0x00, iLen * sizeof(double));
172                     }
173
174                     isDeletable[i] = true;
175                 }
176             }
177
178             for(int i = 0; i < iRows[0] * iCols[0]; i++)
179             {
180                 int iErr    = 0;
181                 double real = 0;
182                 double imag = 0;
183
184                 C2F(wesidu)(pdblInR[0][i], pdblInI[0][i], (piRank[0])+i,
185                             pdblInR[1][i], pdblInI[1][i], (piRank[1])+i,
186                             pdblInR[2][i], pdblInI[2][i], (piRank[2])+i,
187                             &real, &imag, &dblEps, &iErr);
188
189                 if(iErr)
190                 {
191                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "wesidu");
192                     throw iErr;
193                 }
194
195                 pdblOutR[i] = real;
196                 pdblOutI[i] = imag;
197             }
198         }
199     }
200     catch(int error)
201     {
202         iError = error;
203     }
204
205     // free memory
206     for(int i = 0; i < 3; i++)
207     {
208         if(pDblIn[i])
209         {
210             delete pDblIn[i];
211         }
212
213         if(pPoly[i])
214         {
215             delete pPoly[i];
216         }
217
218         if(piRank[i])
219         {
220             delete[] piRank[i];
221         }
222
223         if(pdblInR[i])
224         {
225             delete[] pdblInR[i];
226         }
227
228         if(isDeletable[i])
229         {
230             for(int j = 0; j < iSize; j++)
231             {
232                 delete[] pdblInI[i][j];
233             }
234         }
235
236         if(pdblInI[i])
237         {
238             delete[] pdblInI[i];
239         }
240     }
241
242     /*** retrun output arguments ***/
243     if(iError)
244     {
245         return types::Function::Error;
246     }
247
248     out.push_back(pDblOut);
249     return types::Function::OK;
250 }
251 /*--------------------------------------------------------------------------*/