1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2013 - Scilab Enterprises - Antoine ELIAS
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 //RELEASE 4.0, WGS COPYRIGHT 2000.
17 //
18 //scilab call:
19 // [(x0)(,B(,D))(,V)(,rcnd)] = findBD(jobx0,comuse(,job),A(,B),C(,D),Y
20 // (,U,tol,printw,ldwork))
21 //
22 // [x0,B,V,rcnd] = findBD(1,1,1,A,C,Y,U)
23 // [x0,B,D,V,rcnd] = findBD(1,1,2,A,C,Y,U)
24 // [B,V,rcnd] = findBD(0,1,1,A,C,Y,U)
25 // [B,D,V,rcnd] = findBD(0,1,2,A,C,Y,U)
26 // [x0,V,rcnd] = findBD(1,2,1,A,B,C,Y,U)
27 // [x0,V,rcnd] = findBD(1,2,2,A,B,C,D,Y,U)
28 // [x0,rcnd] = findBD(0,2) (x0 = 0, rcnd = 1)
29 // [x0,V,rcnd] = findBD(1,3,A,C,Y)
30 //
31 //Purpose:
32 // To estimate the initial state and/or the matrices B and D of a
33 // discrete-time system, given the system matrices A, C, and possibly
34 // B, D, and the input and output trajectories of the system.
35 //
36 //Input parameters:
37 // jobx0 - integer option to specify whether or not the initial state
38 // should be computed:
39 // = 1 : compute the initial state x(0);
40 // = 0 : do not compute the initial state (possibly, because
41 // x(0) is known to be zero).
42 // comuse- integer option to specify whether the system matrices B
43 // and D should be computed or used:
44 // = 1 : compute the matrices B and D, as specified by job;
45 // = 2 : use the matrices B and D, as specified by job;
46 // = 3 : do not compute/use the matrices B and D.
47 // job - integer option to determine which of the system matrices
48 // B and D should be computed or used:
49 // = 1 : compute/use the matrix B only (D is known to be zero);
50 // = 2 : compute/use the matrices B and D.
51 // job must not be specified if jobx0 = 0 and comuse = 2, or
52 // if comuse = 3.
53 // A - the n-by-n system state matrix A.
54 // B - the n-by-m system input matrix B (input when jobx0 = 1 and
55 // comuse = 2).
56 // C - the l-by-n system output matrix C.
57 // D - the l-by-m system matrix D (input when jobx0 = 1,
58 // comuse = 2
59 // and job = 2).
60 // Y - the t-by-l output-data sequence matrix. Column j of Y
61 // contains the t values of the j-th output component for
62 // consecutive time increments.
63 // U - the t-by-m input-data sequence matrix (input when jobx0 = 1
64 // and comuse = 2, or comuse = 1). Column j of U
65 // contains the t values of the j-th input component for
66 // consecutive time increments.
67 // tol - (optional) tolerance used for estimating the rank of
68 // matrices. If tol > 0, then the given value of tol is
69 // used as a lower bound for the reciprocal condition number;
70 // an m-by-n matrix whose estimated condition number is less
71 // than 1/tol is considered to be of full rank.
72 // Default: m*n*epsilon_machine where epsilon_machine is the
73 // relative machine precision.
74 // printw- (optional) switch for printing the warning messages.
75 // = 1: print warning messages;
76 // = 0: do not print warning messages.
77 // Default: printw = 0.
78 // ldwork- (optional) the workspace size.
79 // Default : computed by the formula
80 // LDWORK = MAX( minimum workspace size needed, 2*CSIZE/3,
81 // CSIZE - ( m + l )*t - 2*n*( n + m + l ) - l*m )
82 // where CSIZE is the cache size in double precision words.
83 //
84 //Output parameters:
85 // x0 - the n-vector of estimated initial state x(0) (output when
86 // jobx0 = 1, and set to 0 when jobx0 = 0 and comuse = 2).
87 // B - the n-by-m system input matrix B (output when comuse = 1).
88 // D - the l-by-m system matrix D (output when comuse = 1 and
89 // job = 2).
90 // V - the n-by-n orthogonal matrix which reduces A to a real Schur
91 // form (output when jobx0 = 1 or comuse = 1).
92 // rcnd - (optional) vector of length 1 (or 2, if comuse = 1 and
93 // job = 2), containing the reciprocal condition numbers of
94 // the matrices involved in rank decisions.
95 //
96 //Contributor:
97 // V. Sima, Research Institute for Informatics, Bucharest, April 2000.
98 //
99 //Revisions:
100 // V. Sima, July 2000, February 2001.
101 //
102 //**********************************************************************
103 //
104
105 #include <string.h>
106 #include "gw_cacsd.h"
107 #include "api_scilab.h"
108 #include "Scierror.h"
109 #include "sciprint.h"
110 #include "localization.h"
111 #include "BOOL.h"
112 #include "sci_malloc.h"
113 #include "gw_cacsd.h"
114 #include "elem_common.h"
115 /*--------------------------------------------------------------------------*/
116 extern int C2F(ib01cd)(char*, char*, char*, int*, int*, int*, int*, double*, int*, double*,
117 int*, double*, int*, double*, int*, double*, int*, double*, int*, double*, double*,
118 int*, double*, int*, double*, int*, int*, int*);
119 /*--------------------------------------------------------------------------*/
sci_findbd(char * fname,void * pvApiCtx)120 int sci_findbd(char *fname, void* pvApiCtx)
121 {
122 SciErr sciErr;
123 int iOne = 1;
124
125 int iTASK = 0;
126 char cJOBX0 = 'N';
127 int iCUSE = 0;
128 char cCOMUSE = 'N';
129 int iIJOB = 0;
130 char cJOB = 'D';
131
132 int iN = 0;
133 int iM = 0;
134 int iL = 1;
135 int iNSMP = 1;
136 int iNCOL = 0;
137 int iMINSMP = 0;
138 int iIQ = 0;
139
140 int iSizeA = 0;
141 int iSizeB = 0;
142 int iSizeC = 0;
143 int iSizeD = 0;
144 int iSizeY = 0;
145 int iSizeU = 0;
146
147 BOOL bPRINTW = FALSE;
148 double dblTOL = 0;
149
150 double* pdblA = NULL;
151 double* pdblB = NULL;
152 double* pdblC = NULL;
153 double* pdblD = NULL;
154 double* pdblY = NULL;
155 double* pdblU = NULL;
156
157 int iIPS = 0;
158 int iIP = 0;
159 int iRhs = nbInputArgument(pvApiCtx);
160 int iLhs = nbOutputArgument(pvApiCtx);
161 int iCurrentLhs = 1;
162 int iPos = 0;
163
164 int iLDA = 0;
165 int iLDB = 0;
166 int iLDC = 0;
167 int iLDD = 0;
168 int iLDV = 0;
169 int iLDY = 0;
170 int iLDU = 0;
171 int iLIWORK = 0;
172 int iLDWORK = 0;
173 int iLDWMIN = 0;
174
175 int iWARN = 0;
176 int iINFO = 0;
177
178 /*temporary variables*/
179 double* pA = NULL;
180 double* pB = NULL;
181 double* pC = NULL;
182 double* pD = NULL;
183 double* pDWORK = NULL;
184 int* pIWORK = NULL;
185 double* pU = NULL;
186 double* pV = NULL;
187 double* pX0 = NULL;
188 double* pY = NULL;
189
190 CheckInputArgumentAtLeast(pvApiCtx, 2);
191 CheckOutputArgumentAtLeast(pvApiCtx, 1);
192
193
194 // jobx0
195 CHECK_PARAM(pvApiCtx, 1);
196 iTASK = getIntegerValue(pvApiCtx, 1);
197
198 if (iTASK == 1)
199 {
200 cJOBX0 = 'X';
201 }
202
203 // comuse
204 CHECK_PARAM(pvApiCtx, 2);
205 iCUSE = getIntegerValue(pvApiCtx, 2);
206
207 if (iCUSE == 1)
208 {
209 cCOMUSE = 'C';
210 }
211
212 if (iCUSE == 2)
213 {
214 cCOMUSE = 'U';
215 }
216
217 if (iTASK == 1 && iCUSE == 2)
218 {
219 if (iRhs < 5)
220 {
221 Scierror(999, _("%s: Wrong number of input argument(s): At least %d expected.\n"), fname, 5);
222 return 0;
223 }
224 }
225 else if (iCUSE == 1)
226 {
227 if (iRhs < 6)
228 {
229 Scierror(999, _("%s: Wrong number of input argument(s): At least %d expected.\n"), fname, 6);
230 return 0;
231 }
232 }
233
234
235 // job
236 if ((iTASK == 1 && iCUSE == 2) || iCUSE == 1)
237 {
238 CHECK_PARAM(pvApiCtx, 3);
239 iIJOB = getIntegerValue(pvApiCtx, 3);
240 iIP = 4;
241 if (iIJOB == 1)
242 {
243 cJOB = 'B';
244 }
245 }
246 else
247 {
248 cJOB = 'B';
249 iIJOB = 1;
250 iIP = 3;
251 }
252
253 if (iTASK == 1 && iCUSE == 2)
254 {
255 if (iIJOB == 1 && iRhs < 7)
256 {
257 Scierror(999, _("%s: Wrong number of input argument(s): At least %d expected.\n"), fname, 7);
258 return 0;
259 }
260 else if (iIJOB == 2 && iRhs < 8)
261 {
262 Scierror(999, _("%s: Wrong number of input argument(s): At least %d expected.\n"), fname, 8);
263 return 0;
264 }
265 }
266
267 // A(n,n)
268 iIPS = iIP;
269
270 if (iRhs > iIP)
271 {
272 int* piAddr = NULL;
273 int iCols = 0;
274
275 sciErr = getVarAddressFromPosition(pvApiCtx, iIP, &piAddr);
276 if (sciErr.iErr)
277 {
278 printError(&sciErr, 0);
279 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
280 return 0;
281 }
282
283 getMatrixOfDouble(pvApiCtx, piAddr, &iN, &iCols, &pdblA);
284 if (sciErr.iErr)
285 {
286 printError(&sciErr, 0);
287 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
288 return 0;
289 }
290
291 if (iN != iCols)
292 {
293 Scierror(999, _("%s: Wrong size for input argument #%d: A matrix of size %dx%d expected.\n"), fname, iIP, iN, iN);
294 return 0;
295 }
296
297 iIP++;
298 }
299
300 // B(n,m)
301 if (iTASK == 1 && iCUSE == 2 && iRhs >= iIP)
302 {
303 int* piAddr = NULL;
304 int iRows = 0;
305
306 sciErr = getVarAddressFromPosition(pvApiCtx, iIP, &piAddr);
307 if (sciErr.iErr)
308 {
309 printError(&sciErr, 0);
310 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
311 return 0;
312 }
313
314 getMatrixOfDouble(pvApiCtx, piAddr, &iRows, &iM, &pdblB);
315 if (sciErr.iErr)
316 {
317 printError(&sciErr, 0);
318 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
319 return 0;
320 }
321
322 if (iRows != iN)
323 {
324 Scierror(999, _("%s: Wrong size for input argument #%d: A matrix of size %dx%d expected.\n"), fname, iIP, iN, iN);
325 return 0;
326 }
327
328 iIP++;
329 }
330
331 // C(1,n)
332 if (iRhs >= iIP)
333 {
334 int* piAddr = NULL;
335 int iCols = 0;
336
337 sciErr = getVarAddressFromPosition(pvApiCtx, iIP, &piAddr);
338 if (sciErr.iErr)
339 {
340 printError(&sciErr, 0);
341 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
342 return 0;
343 }
344
345 getMatrixOfDouble(pvApiCtx, piAddr, &iL, &iCols, &pdblC);
346 if (sciErr.iErr)
347 {
348 printError(&sciErr, 0);
349 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
350 return 0;
351 }
352
353 if (iCols != iN)
354 {
355 Scierror(999, _("%s: Wrong size for input argument #%d: A matrix with %d columns expected.\n"), fname, iIP, iN);
356 return 0;
357 }
358
359 if (iL <= 0)
360 {
361 Scierror(999, _("%s: The system has no outputs\n"), fname);
362 return 0;
363 }
364
365 iIP++;
366 }
367
368 // D(1,m)
369 if (iTASK == 1 && iCUSE == 2 && iIJOB == 2 && iRhs >= iIP)
370 {
371 int* piAddr = NULL;
372 int iRows = 0;
373 int iCols = 0;
374
375 sciErr = getVarAddressFromPosition(pvApiCtx, iIP, &piAddr);
376 if (sciErr.iErr)
377 {
378 printError(&sciErr, 0);
379 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
380 return 0;
381 }
382
383 getMatrixOfDouble(pvApiCtx, piAddr, &iRows, &iCols, &pdblD);
384 if (sciErr.iErr)
385 {
386 printError(&sciErr, 0);
387 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
388 return 0;
389 }
390
391 if (iRows != iL || iCols != iM)
392 {
393 Scierror(999, _("%s: Wrong size for input argument #%d: A matrix of size %dx%d expected.\n"), fname, iIP, iL, iM);
394 return 0;
395 }
396
397 iIP++;
398 }
399
400 // Y(txp)
401
402 if (iRhs >= iIP)
403 {
404 int* piAddr = NULL;
405 int iCols = 0;
406
407 sciErr = getVarAddressFromPosition(pvApiCtx, iIP, &piAddr);
408 if (sciErr.iErr)
409 {
410 printError(&sciErr, 0);
411 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
412 return 0;
413 }
414
415 getMatrixOfDouble(pvApiCtx, piAddr, &iNSMP, &iCols, &pdblY);
416 if (sciErr.iErr)
417 {
418 printError(&sciErr, 0);
419 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
420 return 0;
421 }
422
423 if (iCols != iL)
424 {
425 Scierror(999, _("%s: Wrong size for input argument #%d: A matrix with %d columns expected.\n"), fname, iIP, iL);
426 return 0;
427 }
428
429 iIP++;
430 }
431
432 // U(txm)
433 if ((iTASK == 1 && iCUSE == 2) || (iCUSE == 1 && iRhs > iIP))
434 {
435 int* piAddr = NULL;
436 int iRows = 0;
437 int iCols = 0;
438
439 sciErr = getVarAddressFromPosition(pvApiCtx, iIP, &piAddr);
440 if (sciErr.iErr)
441 {
442 printError(&sciErr, 0);
443 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
444 return 0;
445 }
446
447 getMatrixOfDouble(pvApiCtx, piAddr, &iRows, &iCols, &pdblU);
448 if (sciErr.iErr)
449 {
450 printError(&sciErr, 0);
451 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
452 return 0;
453 }
454
455 if (iCUSE == 2 && iCols != iM)
456 {
457 Scierror(999, _("%s: Wrong size for input argument #%d: A matrix with %d columns expected.\n"), fname, iIP, iM);
458 return 0;
459 }
460
461 iM = iCols;
462
463 if (iM > 0)
464 {
465 if (iRows != iNSMP)
466 {
467 Scierror(999, _("%s: Incompatible input arguments #%d and #%d: Same row dimensions expected.\n"), fname, iIP - 1, iIP);
468 return 0;
469 }
470 }
471
472 iIP++;
473 }
474
475 if (iCUSE == 1)
476 {
477 iNCOL = iN * iM;
478 if (iTASK == 1)
479 {
480 iNCOL += iN;
481 }
482 iMINSMP = iNCOL;
483
484 if (iIJOB == 2)
485 {
486 iMINSMP += iM;
487 iIQ = iMINSMP;
488 }
489 else if (iTASK == 0)
490 {
491 iIQ = iMINSMP;
492 }
493 else
494 {
495 iIQ = iMINSMP;
496 }
497 }
498 else
499 {
500 iNCOL = iN;
501 if (iTASK == 1)
502 {
503 iMINSMP = iN;
504 }
505 else
506 {
507 iMINSMP = 0;
508 }
509
510 iIQ = iMINSMP;
511 }
512
513
514 if (iNSMP < iMINSMP)
515 {
516 Scierror(999, _("%s: The number of samples should be at least %d"), fname, iMINSMP);
517 return 0;
518 }
519
520 // tol
521
522 if (iRhs >= iIP)
523 {
524 int* piAddr = NULL;
525
526 sciErr = getVarAddressFromPosition(pvApiCtx, iIP, &piAddr);
527 if (sciErr.iErr)
528 {
529 printError(&sciErr, 0);
530 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
531 return 0;
532 }
533
534 getScalarDouble(pvApiCtx, piAddr, &dblTOL);
535 if (sciErr.iErr)
536 {
537 printError(&sciErr, 0);
538 Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, iIP);
539 return 0;
540 }
541
542 iIP++;
543 }
544
545
546 // printw
547 if (iRhs >= iIP)
548 {
549 int iPrintw = 0;
550 CHECK_PARAM(pvApiCtx, iIP);
551 iPrintw = getIntegerValue(pvApiCtx, iIP);
552
553 switch (iPrintw)
554 {
555 case 1 :
556 bPRINTW = TRUE;
557 break;
558 case 0 :
559 bPRINTW = FALSE;
560 break;
561 default :
562 Scierror(999, _("%s: Wrong value for input argument #%d: '%s' or '%s' expected.\n"), fname, 10, "0", "1");
563 return 0;
564 }
565 }
566
567
568 // Determine the lengths of working arrays.
569 // The default value for LDWORK is computed using the formula
570 // LDWORK = MAX( minimum value needed, 2*CSIZE/3,
571 // $ CSIZE - ( M + L )*NSMP - 2*N*( N + M + L ) - L*M )
572 // where CSIZE is the cache size in double precision words.
573 // If LDWORK is specified,
574 // but it is less than the minimum workspace size
575 // needed, that minimum value is used instead.
576
577 iLDA = Max(1, iN);
578 iLDB = iLDA;
579 iLDC = iL;
580 iLDD = iLDC;
581 iLDV = iLDA;
582 iLDY = Max(1, iNSMP);
583 iLDU = 1;
584
585 if (iM > 0 && ((iTASK == 1 && iCUSE == 2) || iCUSE == 1))
586 {
587 iLDU = iLDY;
588 }
589
590 iLIWORK = iNCOL;
591
592 if ((iTASK == 0 && iCUSE != 1) || Max(iN, iM) == 0)
593 {
594 iLDWORK = 2;
595 }
596 else
597 {
598 int iISIZE = 0;
599 int iNCP1 = 0;
600 int iIC = 0;
601 int iMINWLS = 0;
602
603 int iLDW1 = 0;
604 int iLDW2 = 0;
605 int iLDW3 = 0;
606 int iINI = 0;
607
608 if (iIJOB == 2)
609 {
610 iLIWORK = Max(iLIWORK, iM);
611 }
612
613 iIQ *= iL;
614 iNCP1 = iNCOL + 1;
615 iISIZE = iNSMP * iL * iNCP1;
616 if (iN > 0 && iTASK == 1)
617 {
618 if (iCUSE == 1)
619 {
620 iIC = 2 * iN * iN + iN;
621 }
622 else
623 {
624 iIC = 2 * iN * iN;
625 }
626 }
627 else
628 {
629 iIC = 0;
630 }
631
632 iMINSMP = iNCOL * iNCP1;
633
634 if (iCUSE == 1)
635 {
636 int iIA = 0;
637
638 if (iIJOB == 2)
639 {
640 iMINWLS += iL * iM * iNCP1;
641 }
642
643 if (iM > 0 && iIJOB == 2)
644 {
645 iIA = iM + Max(2 * iNCOL, iM);
646 }
647 else
648 {
649 iIA = 2 * iNCOL;
650 }
651
652 iLDW1 = iN * iN * iM + Max(iIC, iIA);
653 if (iTASK == 1)
654 {
655 iLDW1 += iL * iN;
656 }
657
658 iLDW2 = iISIZE + Max(iN + Max(iIC, iIA), 6 * iNCOL);
659 iLDW3 = iMINWLS + Max(iIQ * iNCP1 + iLDW1, 6 * iNCOL);
660 if (iM > 0 && iIJOB == 2)
661 {
662 iLDW2 = Max(iLDW2, iISIZE + 2 * iM * iM + 6 * iM);
663 iLDW3 = Max(iLDW3, iMINWLS + 2 * iM * iM + 6 * iM);
664 iINI = 3;
665 }
666 else
667 {
668 iINI = 2;
669 }
670 }
671 else
672 {
673 int iITAU = iIC + iL * iN;
674 iLDW2 = iISIZE + 2 * iN + Max(iIC, 4 * iN);
675 iLDW3 = iMINWLS + 2 * iN + Max(iIQ + iNCP1 + iITAU, 4 * iN);
676 iINI = 2;
677 }
678
679 iLDWMIN = iINI + iN * (iN + iM + iL) + Max(5 * iN, Max(iINI, Min(iLDW2, iLDW3)));
680 iLDWORK = Max(Max(iLDWMIN, 2 * CSIZE / 3), CSIZE - (iM + iL ) * iNSMP - 2 * iN * (iN + iM + iL ) - iL * iM);
681 }
682
683 // ldwork
684 if (iRhs > iIP)
685 {
686 CHECK_PARAM(pvApiCtx, iIP);
687 iLDWORK = Max(getIntegerValue(pvApiCtx, iIP), iLDWMIN);
688 }
689
690 // Allocate variable dimension local arrays.
691 iSizeA = iLDA * iN;
692 iSizeB = iLDB * iM;
693 iSizeC = iLDC * iN;
694 iSizeD = iLDD * iM;
695 iSizeU = iLDU * iM;
696 iSizeY = iLDY * iL;
697
698 pA = (double*)MALLOC(iLDA * iN * sizeof(double));
699 pB = (double*)MALLOC(iLDB * iM * sizeof(double));
700 pC = (double*)MALLOC(iLDC * iN * sizeof(double));
701 pD = (double*)MALLOC(iLDD * iM * sizeof(double));
702 pDWORK = (double*)MALLOC(iLDWORK * sizeof(double));
703 pIWORK = (int*)MALLOC(iLIWORK * sizeof(int));
704 pU = (double*)MALLOC(iLDU * iM * sizeof(double));
705 pV = (double*)MALLOC(iLDV * iN * sizeof(double));
706 pX0 = (double*)MALLOC(iN * sizeof(double));
707 pY = (double*)MALLOC(iLDY * iL * sizeof(double));
708
709 //Copy inputs from scilab workspace to locally allocated arrays.
710
711 iIP = iIPS;
712 if (iRhs >= iIP)
713 {
714 C2F(dcopy)(&iSizeA, pdblA, &iOne, pA, &iOne);
715 if (iTASK == 1 && iCUSE == 2)
716 {
717 if (iN > 0)
718 {
719 C2F(dcopy)(&iSizeB, pdblB, &iOne, pB, &iOne);
720 }
721 }
722
723 C2F(dcopy)(&iSizeC, pdblC, &iOne, pC, &iOne);
724 if (iTASK == 1 && iCUSE == 2 && iIJOB == 2)
725 {
726 C2F(dcopy)(&iSizeD, pdblD, &iOne, pD, &iOne);
727 }
728
729 C2F(dcopy)(&iSizeY, pdblY, &iOne, pY, &iOne);
730 if (iCUSE == 1 || (iTASK == 1 && iCUSE == 2))
731 {
732 C2F(dcopy)(&iSizeU, pdblU, &iOne, pU, &iOne);
733 }
734 }
735
736 //Do the actual computations.
737 C2F(ib01cd)(&cJOBX0 , &cCOMUSE, &cJOB, &iN, &iM, &iL, &iNSMP, pA, &iLDA, pB, &iLDB,
738 pC, &iLDC , pD, &iLDD, pU, &iLDU, pY, &iLDY, pX0, pV, &iLDV, &dblTOL, pIWORK,
739 pDWORK, &iLDWORK, &iWARN, &iINFO);
740
741 if (iWARN != 0 && bPRINTW)
742 {
743 sciprint("IWARN = %d on exit from IB01CD\n", iWARN);
744 }
745
746 if (iINFO != 0)
747 {
748 Scierror(999, _("%s: INFO = %d on exit from IB01CD\n"), fname, iINFO);
749 FREE(pA);
750 FREE(pB);
751 FREE(pC);
752 FREE(pD);
753 FREE(pDWORK);
754 FREE(pIWORK);
755 FREE(pU);
756 FREE(pV);
757 FREE(pX0);
758 FREE(pY);
759 return 0;
760 }
761
762 iCurrentLhs = iRhs + 1;
763 iPos = 1;
764 //Copy output to scilab workspace.
765 if (iTASK == 1 || (iTASK == 0 && iCUSE == 2))
766 {
767 createMatrixOfDouble(pvApiCtx, iCurrentLhs, iN, 1, pX0);
768 AssignOutputVariable(pvApiCtx, iPos) = iCurrentLhs;
769 iCurrentLhs++;
770 iPos++;
771 }
772
773 if (iCUSE == 1 && iLhs >= iPos)
774 {
775 createMatrixOfDouble(pvApiCtx, iCurrentLhs, iN, iM, pB);
776 AssignOutputVariable(pvApiCtx, iPos) = iCurrentLhs;
777 iCurrentLhs++;
778 iPos++;
779
780 if (iLhs >= iPos)
781 {
782 createMatrixOfDouble(pvApiCtx, iCurrentLhs, iL, iM, pD);
783 AssignOutputVariable(pvApiCtx, iPos) = iCurrentLhs;
784 iCurrentLhs++;
785 iPos++;
786 }
787 }
788
789 if ((iTASK == 1 || iCUSE == 1) && iLhs >= iPos)
790 {
791 createMatrixOfDouble(pvApiCtx, iCurrentLhs, iN, iN, pV);
792 AssignOutputVariable(pvApiCtx, iPos) = iCurrentLhs;
793 iCurrentLhs++;
794 iPos++;
795 }
796
797 if (iLhs >= iPos)
798 {
799 int iNO = 1;
800 if (iCUSE == 1 && iM > 0 && iIJOB == 2)
801 {
802 iNO = 2;
803 }
804
805 createMatrixOfDouble(pvApiCtx, iCurrentLhs, iNO, 1, pDWORK + 1);
806 AssignOutputVariable(pvApiCtx, iPos) = iCurrentLhs;
807 }
808
809 ReturnArguments(pvApiCtx);
810
811 FREE(pA);
812 FREE(pB);
813 FREE(pC);
814 FREE(pD);
815 FREE(pDWORK);
816 FREE(pIWORK);
817 FREE(pU);
818 FREE(pV);
819 FREE(pX0);
820 FREE(pY);
821 return 0;
822 }
823