1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2012 - 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 #include <string.h>
16 #include "gw_cacsd.h"
17 #include "sci_rankqr.h"
18 #include "api_scilab.h"
19 #include "Scierror.h"
20 #include "localization.h"
21
22 extern int C2F(riccsl)();
23 extern int C2F(riccms)();
24 extern int C2F(ricdsl)();
25 extern int C2F(ricdmf)();
26
sci_ricc(char * fname,void * pvApiCtx)27 int sci_ricc(char *fname, void* pvApiCtx)
28 {
29 SciErr sciErr;
30
31 int* piAddrlA = NULL;
32 double* lA = NULL;
33 int* piAddrlD = NULL;
34 double* lD = NULL;
35 int* piAddrlC = NULL;
36 double* lC = NULL;
37 char* lTYPE = NULL;
38 char* lMETHOD = NULL;
39 double* lX = NULL;
40 double* lWR = NULL;
41 double* lWI = NULL;
42 double* lRCOND = NULL;
43 double* lFERR = NULL;
44 int* lIWORK = NULL;
45 int* lBWORK = NULL;
46 double* lDWORK = NULL;
47
48 int* piAddrlTYPE = NULL;
49 int* piAddrlMETHOD = NULL;
50
51 BOOL WANTC = 0, WANTD = 0, WSCHUR = 0, WSIGN = 0, WINVF = 0;
52
53 int minrhs = 4;
54 int maxrhs = 5;
55 int minLhs = 0;
56 int maxLhs = 3;
57
58 int N = 0, LWORKMIN = 0, INFO = 0;
59 int MA = 0, NA = 0;
60 int MC = 0, NC = 0;
61 int MD = 0, ND = 0;
62 int M1 = 0, N1 = 0;
63
64 int iOne = 1, iSize = 0, k = 0;
65
66 CheckInputArgument(pvApiCtx, minrhs, maxrhs);
67 CheckOutputArgument(pvApiCtx, minLhs, maxLhs);
68
69 sciErr = getVarAddressFromPosition(pvApiCtx, 1, &piAddrlA);
70 if (sciErr.iErr)
71 {
72 printError(&sciErr, 0);
73 return 1;
74 }
75
76 // Retrieve a matrix of double at position 1.
77 sciErr = getMatrixOfDouble(pvApiCtx, piAddrlA, &MA, &NA, &lA);
78 if (sciErr.iErr)
79 {
80 printError(&sciErr, 0);
81 Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 1);
82 return 1;
83 }
84
85 if (MA != NA)
86 {
87 Scierror(999, _("%s: A must be a square matrix.\n"), fname);
88 return 1;
89 }
90
91 sciErr = getVarAddressFromPosition(pvApiCtx, 2, &piAddrlD);
92 if (sciErr.iErr)
93 {
94 printError(&sciErr, 0);
95 return 1;
96 }
97
98 // Retrieve a matrix of double at position 2.
99 sciErr = getMatrixOfDouble(pvApiCtx, piAddrlD, &MD, &ND, &lD);
100 if (sciErr.iErr)
101 {
102 printError(&sciErr, 0);
103 Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 2);
104 return 1;
105 }
106
107 if (MD != ND)
108 {
109 Scierror(999, _("%s: D must be a square matrix.\n"), fname);
110 return 1;
111 }
112
113 sciErr = getVarAddressFromPosition(pvApiCtx, 3, &piAddrlC);
114 if (sciErr.iErr)
115 {
116 printError(&sciErr, 0);
117 return 1;
118 }
119
120 // Retrieve a matrix of double at position 3.
121 sciErr = getMatrixOfDouble(pvApiCtx, piAddrlC, &MC, &NC, &lC);
122 if (sciErr.iErr)
123 {
124 printError(&sciErr, 0);
125 Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 3);
126 return 1;
127 }
128
129 if (MC != NC)
130 {
131 Scierror(999, _("%s: C must be a square matrix.\n"), fname);
132 return 1;
133 }
134
135 if (MA != MC || MC != MD || MA != MD)
136 {
137 Scierror(999, _("%s: The matrices A, C and D must have the same order.\n"), fname);
138 return 1;
139 }
140
141 N = MA;
142
143 sciErr = getVarAddressFromPosition(pvApiCtx, 4, &piAddrlTYPE);
144 if (sciErr.iErr)
145 {
146 printError(&sciErr, 0);
147 return 1;
148 }
149
150 // Retrieve a matrix of double at position 4.
151 if (getAllocatedSingleString(pvApiCtx, piAddrlTYPE, &lTYPE))
152 {
153 Scierror(202, _("%s: Wrong type for argument #%d: string expected.\n"), fname, 4);
154 return 1;
155 }
156
157 WANTC = (strcmp((lTYPE), "cont") == 0 || strcmp((lTYPE), "CONT") == 0);
158 WANTD = (strcmp((lTYPE), "disc") == 0 || strcmp((lTYPE), "DISC") == 0);
159 freeAllocatedSingleString(lTYPE);
160
161 if (WANTC == FALSE && WANTD == FALSE)
162 {
163 Scierror(999, _("%s: Wrong value for input argument #%d: Type must be continuous or discrete.\n"), fname, 4);
164 return 1;
165 }
166
167 k = 5;
168 WSCHUR = TRUE;
169 if (nbInputArgument(pvApiCtx) == 5)
170 {
171 sciErr = getVarAddressFromPosition(pvApiCtx, 5, &piAddrlMETHOD);
172 if (sciErr.iErr)
173 {
174 printError(&sciErr, 0);
175 return 1;
176 }
177
178 // Retrieve a matrix of double at position 5.
179 if (getAllocatedSingleString(pvApiCtx, piAddrlMETHOD, &lMETHOD))
180 {
181 Scierror(202, _("%s: Wrong type for argument #%d: string expected.\n"), fname, 5);
182 return 1;
183 }
184
185 WSCHUR = (strcmp((lMETHOD), "schr") == 0 || strcmp((lMETHOD), "SCHR") == 0);
186 if (WANTC)
187 {
188 WSIGN = (strcmp((lMETHOD), "sign") == 0 || strcmp((lMETHOD), "SIGN") == 0);
189 if (WSCHUR == FALSE && WSIGN == FALSE)
190 {
191 freeAllocatedSingleString(lMETHOD);
192 Scierror(999, _("%s: Wrong value for input argument #%d: Method must be schur or sign.\n"), fname, 5);
193 return 1;
194 }
195 }
196 else
197 {
198 WINVF = (strcmp((lMETHOD), "invf") == 0 || strcmp((lMETHOD), "INVF") == 0);
199 if (WSCHUR == FALSE && WINVF == FALSE)
200 {
201 freeAllocatedSingleString(lMETHOD);
202 Scierror(999, _("%s: Wrong value for input argument #%d: Method must be schur or invf.\n"), fname, 5);
203 return 1;
204 }
205 }
206
207 k = 6;
208 freeAllocatedSingleString(lMETHOD);
209 }
210
211 sciErr = allocMatrixOfDouble(pvApiCtx, k, N, N, &lX);
212 if (sciErr.iErr)
213 {
214 printError(&sciErr, 0);
215 Scierror(999, _("%s: Memory allocation error.\n"), fname);
216 return 1;
217 }
218
219 sciErr = allocMatrixOfDouble(pvApiCtx, k + 1, N, iOne, &lWR);
220 if (sciErr.iErr)
221 {
222 printError(&sciErr, 0);
223 Scierror(999, _("%s: Memory allocation error.\n"), fname);
224 return 1;
225 }
226
227 sciErr = allocMatrixOfDouble(pvApiCtx, k + 2, N, iOne, &lWI);
228 if (sciErr.iErr)
229 {
230 printError(&sciErr, 0);
231 Scierror(999, _("%s: Memory allocation error.\n"), fname);
232 return 1;
233 }
234
235 sciErr = allocMatrixOfDouble(pvApiCtx, k + 3, iOne, iOne, &lRCOND);
236 if (sciErr.iErr)
237 {
238 printError(&sciErr, 0);
239 Scierror(999, _("%s: Memory allocation error.\n"), fname);
240 return 1;
241 }
242
243 sciErr = allocMatrixOfDouble(pvApiCtx, k + 4, iOne, iOne, &lFERR);
244 if (sciErr.iErr)
245 {
246 printError(&sciErr, 0);
247 Scierror(999, _("%s: Memory allocation error.\n"), fname);
248 return 1;
249 }
250
251 iSize = Max(2 * N, N * N);
252 sciErr = allocMatrixOfDoubleAsInteger(pvApiCtx, k + 5, iOne, iSize, &lIWORK);
253 if (sciErr.iErr)
254 {
255 printError(&sciErr, 0);
256 Scierror(999, _("%s: Memory allocation error.\n"), fname);
257 return 1;
258 }
259
260 iSize = 2 * N;
261 sciErr = allocMatrixOfDoubleAsInteger(pvApiCtx, k + 6, iOne, iSize, &lBWORK);
262 if (sciErr.iErr)
263 {
264 printError(&sciErr, 0);
265 Scierror(999, _("%s: Memory allocation error.\n"), fname);
266 return 1;
267 }
268
269 if (WANTC)
270 {
271 if (WSCHUR)
272 {
273 LWORKMIN = 9 * N * N + 4 * N + Max(1, 6 * N);
274 }
275 else if (WSIGN)
276 {
277 LWORKMIN = 9 * N * N + 7 * N + 1;
278 }
279 }
280 else
281 {
282 if (WSCHUR)
283 {
284 LWORKMIN = 12 * N * N + 22 * N + Max(16, 4 * N);
285 }
286 else if (WINVF)
287 {
288 LWORKMIN = 28 * N * N + 2 * N + Max(1, 2 * N);
289 }
290 }
291
292 sciErr = allocMatrixOfDouble(pvApiCtx, k + 7, iOne, LWORKMIN, &lDWORK);
293 if (sciErr.iErr)
294 {
295 printError(&sciErr, 0);
296 Scierror(999, _("%s: Memory allocation error.\n"), fname);
297 return 1;
298 }
299
300 if (WANTC)
301 {
302 if (WSCHUR)
303 {
304 C2F(riccsl)("N", &N, lA, &N, "U", lC, &N, lD,
305 &N, lX, &N, lWR, lWI, lRCOND,
306 lFERR, lDWORK, &LWORKMIN, lIWORK,
307 lBWORK, &INFO);
308
309 if (INFO != 0)
310 {
311 Scierror(999, _("%s: RICCSL exit with info = %d.\n"), fname, INFO);
312 return 1;
313 }
314 }
315 else if (WSIGN)
316 {
317 C2F(riccms)("N", &N, lA, &N, "U", lC, &N, lD,
318 &N, lX, &N, lWR, lWI, lRCOND,
319 lFERR, lDWORK, &LWORKMIN, lIWORK, &INFO);
320
321 if (INFO != 0)
322 {
323 Scierror(999, _("%s: RICCMS exit with info = %d.\n"), fname, INFO);
324 return 1;
325 }
326 }
327 }
328 else
329 {
330 if (WSCHUR)
331 {
332 C2F(ricdsl)("N", &N, lA, &N, "U", lC, &N, lD,
333 &N, lX, &N, lWR, lWI, lRCOND,
334 lFERR, lDWORK, &LWORKMIN, lIWORK,
335 lBWORK, &INFO);
336
337 if (INFO != 0)
338 {
339 Scierror(999, _("%s: RICDSL exit with info = %d.\n"), fname, INFO);
340 return 1;
341 }
342 }
343 else if (WINVF)
344 {
345 C2F(ricdmf)("N", &N, lA, &N, "U", lC, &N, lD,
346 &N, lX, &N, lWR, lWI, lRCOND,
347 lFERR, lDWORK, &LWORKMIN, lIWORK, &INFO);
348
349 if (INFO != 0)
350 {
351 Scierror(999, _("%s: RICDMF exit with info = %d.\n"), fname, INFO);
352 return 1;
353 }
354 }
355 }
356
357 if (nbOutputArgument(pvApiCtx) <= 1)
358 {
359 AssignOutputVariable(pvApiCtx, 1) = k;
360 }
361 else if (nbOutputArgument(pvApiCtx) == 2)
362 {
363 AssignOutputVariable(pvApiCtx, 1) = k;
364 AssignOutputVariable(pvApiCtx, 2) = k + 3;
365 }
366 else if (nbOutputArgument(pvApiCtx) == 3)
367 {
368 AssignOutputVariable(pvApiCtx, 1) = k;
369 AssignOutputVariable(pvApiCtx, 2) = k + 3;
370 AssignOutputVariable(pvApiCtx, 3) = k + 4;
371 }
372
373 ReturnArguments(pvApiCtx);
374 return 0;
375 }
376