ad38e50670af5840757fa5cbcbff1b1db7b1cea8
[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     bool isDeletable = false;
47
48     if (in.size() != 3)
49     {
50         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ppol", 3);
51         return types::Function::Error;
52     }
53
54     if (_iRetCount > 1)
55     {
56         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "ppol", 1);
57         return types::Function::Error;
58     }
59
60     /*** get inputs arguments ***/
61     // get poles
62     if (in[2]->isDouble() == false)
63     {
64         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 3);
65         return types::Function::Error;
66     }
67
68     pDblP = in[2]->getAs<types::Double>();
69     iSizeP = pDblP->getSize();
70
71     // get matrix B
72     if (in[1]->isDouble() == false)
73     {
74         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 2);
75         return types::Function::Error;
76     }
77
78     pDblB = in[1]->clone()->getAs<types::Double>();
79
80     if (pDblB->isComplex())
81     {
82         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ppol", 2);
83         return types::Function::Error;
84     }
85
86     iColB = pDblB->getCols();
87
88     // get matrix A
89     if (in[0]->isDouble() == false)
90     {
91         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 1);
92         return types::Function::Error;
93     }
94
95     pDblA = in[0]->clone()->getAs<types::Double>();
96
97     if (pDblA->isComplex())
98     {
99         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ppol", 1);
100         return types::Function::Error;
101     }
102
103     if (pDblA->getRows() != pDblA->getCols())
104     {
105         Scierror(999, _("%s: Wrong size for input argument #%d: A square matrix expected.\n"), "ppol", 1);
106         return types::Function::Error;
107     }
108
109     if (pDblA->getRows() != iSizeP || pDblB->getRows() != iSizeP)
110     {
111         Scierror(999, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "ppol");
112         return types::Function::Error;
113     }
114
115     if (iSizeP == 0)
116     {
117         out.push_back(types::Double::Empty());
118         return types::Function::OK;
119     }
120
121     /*** perform operations ***/
122     types::Double* pDblOut = new types::Double(iColB, iSizeP);
123
124     double* pdblG   = pDblOut->get();
125     double* pdblZ   = new double[iSizeP * iSizeP];
126     int iSizeWork   = std::max(iSizeP * iColB + 3 * iColB,
127                                iColB * iColB + iColB * std::max(2, iColB) + 3 * iColB + 2 * iSizeP);
128     double* pdblW   = new double[iSizeWork];
129     int* piW        = new int[iSizeP];
130
131     double* pdblPReal = pDblP->get();
132     double* pdblPImg  = NULL;
133     if (pDblP->isComplex())
134     {
135         pdblPImg = pDblP->getImg();
136     }
137     else
138     {
139         pdblPImg = new double[iSizeP];
140         memset(pdblPImg, 0x00, iSizeP * sizeof(double));
141         isDeletable = true;
142     }
143
144     int idc = 0;
145     int inc = 0;
146     // calcul de la forme canonique orthogonale
147     C2F(ssxmc)(&iSizeP, &iColB, pDblA->get(), &iSizeP, pDblB->get(), &idc, &inc, piW, pdblZ,
148                pdblW + iColB, pdblW + iColB + iSizeP * iColB, pdblW + 2 * iColB + iSizeP * iColB,
149                pdblW, &dblEps, &iOne);
150
151     if (idc < iSizeP)
152     {
153         Scierror(999, _("%s: Uncontrollable system.\n"), "ppol");
154         delete[] pdblZ;
155         delete[] pdblW;
156         delete[] piW;
157         if (isDeletable)
158         {
159             delete[] pdblPImg;
160         }
161         delete pDblA;
162         delete pDblB;
163         pDblOut->killMe();
164         return types::Function::Error;
165     }
166
167     double* pdblW1 = pdblW  + iColB + iColB * iColB;
168     double* pdblW2 = pdblW1 + iColB * std::max(2, iColB);
169     double* pdblW3 = pdblW2 + iSizeP;
170     double* pdblW4 = pdblW3 + iSizeP;
171     double* pdblW5 = pdblW4 + iColB;
172
173     // placement de poles
174     int iErr = 0;
175     C2F(polmc)( &iSizeP, &iColB, &iSizeP, &iColB, pDblA->get(), pDblB->get(), pdblG, pdblPReal, pdblPImg,
176                 pdblZ, &inc, piW, &iErr, pdblW, pdblW + iColB, pdblW1, pdblW2, pdblW3, pdblW4, pdblW5);
177
178     if (iErr)
179     {
180         Scierror(999, _("%s: Uncontrollable system.\n"), "ppol");
181         delete[] pdblZ;
182         delete[] pdblW;
183         delete[] piW;
184         if (isDeletable)
185         {
186             delete[] pdblPImg;
187         }
188         delete pDblA;
189         delete pDblB;
190         pDblOut->killMe();
191         return types::Function::Error;
192     }
193
194     // free memory
195     delete[] pdblZ;
196     delete[] pdblW;
197     delete[] piW;
198     if (isDeletable)
199     {
200         delete[] pdblPImg;
201     }
202     delete pDblA;
203     delete pDblB;
204
205     /*** retrun output arguments ***/
206     out.push_back(pDblOut);
207     return types::Function::OK;
208 }
209 /*--------------------------------------------------------------------------*/