* Bug 16160 fixed: ppol changed values of third input variable
[scilab.git] / scilab / modules / cacsd / sci_gateway / cpp / sci_ppol.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 extern "C"
20 {
21 #include "Scierror.h"
22 #include "localization.h"
23 #include "elem_common.h"
24
25     extern void C2F(ssxmc)( int*, int*, double*, int*, double*, int*, int*, int*, double*,
26                             double*, double*, double*, double*, double*, int*);
27
28     extern void C2F(polmc)( int*, int*, int*, int*, double*, double*, double*, double*, double*, double*,
29                             int*, int*, int*, double*, double*, double*, double*, double*, double*, double*);
30
31 }
32
33 /*--------------------------------------------------------------------------*/
34 types::Function::ReturnValue sci_ppol(types::typed_list &in, int _iRetCount, types::typed_list &out)
35 {
36     types::Double* pDblA = NULL;
37     types::Double* pDblB = NULL;
38     types::Double* pDblP = NULL;
39
40     double dblEps = 0.1 * sqrt(nc_eps_machine());
41
42     int iOne   = 1;
43     int iSizeP = 0;
44     int iColB  = 0;
45
46     if (in.size() != 3)
47     {
48         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ppol", 3);
49         return types::Function::Error;
50     }
51
52     if (_iRetCount > 1)
53     {
54         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "ppol", 1);
55         return types::Function::Error;
56     }
57
58     /*** get inputs arguments ***/
59     // get poles
60     if (in[2]->isDouble() == false)
61     {
62         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 3);
63         return types::Function::Error;
64     }
65
66     pDblP = in[2]->getAs<types::Double>();
67     iSizeP = pDblP->getSize();
68
69     // get matrix B
70     if (in[1]->isDouble() == false)
71     {
72         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 2);
73         return types::Function::Error;
74     }
75
76     pDblB = in[1]->getAs<types::Double>();
77
78     if (pDblB->isComplex())
79     {
80         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ppol", 2);
81         return types::Function::Error;
82     }
83
84     iColB = pDblB->getCols();
85
86     // get matrix A
87     if (in[0]->isDouble() == false)
88     {
89         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 1);
90         return types::Function::Error;
91     }
92
93     pDblA = in[0]->getAs<types::Double>();
94
95     if (pDblA->isComplex())
96     {
97         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ppol", 1);
98         return types::Function::Error;
99     }
100
101     if (pDblA->getRows() != pDblA->getCols())
102     {
103         Scierror(999, _("%s: Wrong size for input argument #%d: A square matrix expected.\n"), "ppol", 1);
104         return types::Function::Error;
105     }
106
107     if (pDblA->getRows() != iSizeP || pDblB->getRows() != iSizeP)
108     {
109         Scierror(999, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "ppol");
110         return types::Function::Error;
111     }
112
113     if (iSizeP == 0)
114     {
115         out.push_back(types::Double::Empty());
116         return types::Function::OK;
117     }
118
119     /*** perform operations ***/
120
121     pDblA = pDblA->clone();
122     pDblB = pDblB->clone();
123     pDblP = pDblP->clone();
124     pDblP->setComplex(true);
125
126     types::Double* pDblOut = new types::Double(iColB, iSizeP);
127
128     double* pdblG   = pDblOut->get();
129     double* pdblZ   = new double[iSizeP * iSizeP];
130     int iSizeWork   = std::max(iSizeP * iColB + 3 * iColB,
131                                iColB * iColB + iColB * std::max(2, iColB) + 3 * iColB + 2 * iSizeP);
132     double* pdblW   = new double[iSizeWork];
133     int* piW        = new int[iSizeP];
134
135     double* pdblPReal = pDblP->get();
136     double* pdblPImg  = pDblP->getImg();
137
138     int idc = 0;
139     int inc = 0;
140     // calcul de la forme canonique orthogonale
141     C2F(ssxmc)(&iSizeP, &iColB, pDblA->get(), &iSizeP, pDblB->get(), &idc, &inc, piW, pdblZ,
142                pdblW + iColB, pdblW + iColB + iSizeP * iColB, pdblW + 2 * iColB + iSizeP * iColB,
143                pdblW, &dblEps, &iOne);
144
145     if (idc < iSizeP)
146     {
147         Scierror(999, _("%s: Uncontrollable system.\n"), "ppol");
148         delete[] pdblZ;
149         delete[] pdblW;
150         delete[] piW;
151         pDblA->killMe();
152         pDblB->killMe();
153         pDblP->killMe();
154         pDblOut->killMe();
155         return types::Function::Error;
156     }
157
158     double* pdblW1 = pdblW  + iColB + iColB * iColB;
159     double* pdblW2 = pdblW1 + iColB * std::max(2, iColB);
160     double* pdblW3 = pdblW2 + iSizeP;
161     double* pdblW4 = pdblW3 + iSizeP;
162     double* pdblW5 = pdblW4 + iColB;
163
164     // placement de poles
165     int iErr = 0;
166     C2F(polmc)( &iSizeP, &iColB, &iSizeP, &iColB, pDblA->get(), pDblB->get(), pdblG, pdblPReal, pdblPImg,
167                 pdblZ, &inc, piW, &iErr, pdblW, pdblW + iColB, pdblW1, pdblW2, pdblW3, pdblW4, pdblW5);
168
169     // free memory
170     delete[] pdblZ;
171     delete[] pdblW;
172     delete[] piW;
173     pDblA->killMe();
174     pDblB->killMe();
175     pDblP->killMe();
176
177     if (iErr)
178     {
179         Scierror(999, _("%s: Uncontrollable system.\n"), "ppol");
180         pDblOut->killMe();
181         return types::Function::Error;
182     }
183
184     /*** return output arguments ***/
185     out.push_back(pDblOut);
186     return types::Function::OK;
187 }
188 /*--------------------------------------------------------------------------*/