* Bug #7794 fixed - Example of findABCD help page failed.
[scilab.git] / scilab / modules / cacsd / sci_gateway / c / sci_findbd.c
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 /*--------------------------------------------------------------------------*/
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, 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, 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, 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, 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, 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, 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 }