7e31a43d661032dd6988be2c2c2c6c99707bf9f9
[scilab.git] / scilab / modules / optimization / sci_gateway / cpp / sci_semidef.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2013 - 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 "optimization_gw.hxx"
17 #include "function.hxx"
18 #include "double.hxx"
19
20 extern "C"
21 {
22 #include "localization.h"
23 #include "Scierror.h"
24     extern void C2F(spf)(int*, int*, double*, int*, double*, double*, double*,
25                          double*, double*, double*, double*, double*, int*,
26                          double*, int*, int*, int*);
27 }
28
29
30 /*--------------------------------------------------------------------------*/
31 types::Function::ReturnValue sci_semidef(types::typed_list &in, int _iRetCount, types::typed_list &out)
32 {
33     types::Double* pDblX    = NULL;
34     types::Double* pDblZ    = NULL;
35     types::Double* pDblF    = NULL;
36     types::Double* pDblB    = NULL;
37     types::Double* pDblC    = NULL;
38     types::Double* pDblOpt  = NULL;
39
40     int* piB = NULL;
41
42     double dNu      = 0;
43     double dAbstol  = 0;
44     double dReltol  = 0;
45     double dTv      = 0;
46     double dIter    = 0;
47
48     int iSizeX = 0;
49     int iSizeB = 0;
50     int iIter  = 0;
51     int iInfo  = 0;
52
53     if (in.size() < 1 || in.size() > 6)
54     {
55         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "semidef", 1, 6);
56         return types::Function::Error;
57     }
58
59     if (_iRetCount > 4)
60     {
61         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "semidef", 1, 4);
62         return types::Function::Error;
63     }
64
65     /*** get inputs arguments ***/
66     // get X
67     if (in[0]->isDouble() == false)
68     {
69         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 1);
70         return types::Function::Error;
71     }
72
73     pDblX = in[0]->clone()->getAs<types::Double>();
74
75     if (pDblX->isComplex())
76     {
77         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 1);
78         return types::Function::Error;
79     }
80
81     iSizeX = pDblX->getSize();
82
83     // get Z
84     if (in[1]->isDouble() == false)
85     {
86         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 2);
87         return types::Function::Error;
88     }
89
90     pDblZ = in[1]->clone()->getAs<types::Double>();
91
92     if (pDblZ->isComplex())
93     {
94         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 2);
95         return types::Function::Error;
96     }
97
98     // get F
99     if (in[2]->isDouble() == false)
100     {
101         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 3);
102         return types::Function::Error;
103     }
104
105     pDblF = in[2]->getAs<types::Double>();
106
107     if (pDblF->isComplex())
108     {
109         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 3);
110         return types::Function::Error;
111     }
112
113     // get blocksize
114     if (in[3]->isDouble() == false)
115     {
116         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 4);
117         return types::Function::Error;
118     }
119
120     pDblB = in[3]->getAs<types::Double>();
121
122     if (pDblB->isComplex())
123     {
124         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 4);
125         return types::Function::Error;
126     }
127
128     double* pdblB = pDblB->get();
129     iSizeB = pDblB->getSize();
130     piB = new int[iSizeB];
131     for (int i = 0; i < iSizeB; i++)
132     {
133         piB[i] = (int)pdblB[i];
134     }
135
136     // get C
137     if (in[4]->isDouble() == false)
138     {
139         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 5);
140         return types::Function::Error;
141     }
142
143     pDblC = in[4]->getAs<types::Double>();
144
145     if (pDblC->isComplex())
146     {
147         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 5);
148         return types::Function::Error;
149     }
150
151     if (pDblC->getRows()* pDblC->getCols() != pDblX->getRows()*pDblX->getCols())
152     {
153         Scierror(999, _("%s: Wrong size for input argument #%d: Input argument %d size expected.\n"), "semidef", 5, 1);
154         return types::Function::Error;
155     }
156
157     // get Option
158     if (in[5]->isDouble() == false)
159     {
160         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 6);
161         return types::Function::Error;
162     }
163
164     pDblOpt = in[5]->getAs<types::Double>();
165
166     if (pDblOpt->isComplex())
167     {
168         Scierror(999, _("%s: Wrong type for input argument #%d: A real vector expected.\n"), "semidef", 6);
169         return types::Function::Error;
170     }
171
172     if (pDblOpt->getRows() != 1 || pDblOpt->getCols() != 5)
173     {
174         Scierror(999, _("%s: Wrong size for input argument #%d: A vector of size %d x %d expected.\n"), "semidef", 6, 1, 5);
175         return types::Function::Error;
176     }
177
178     dNu      = pDblOpt->get(0);
179     dAbstol  = pDblOpt->get(1);
180     dReltol  = pDblOpt->get(2);
181     dTv      = pDblOpt->get(3);
182     iIter    = (int)pDblOpt->get(4);
183
184     /*** perform operations ***/
185     types::Double* pDblUl = new types::Double(1, 2);
186     int* piWork = new int[iSizeX];
187
188     // compute size of working table
189     int n       = 0;
190     int sz      = 0;
191     int upsz    = 0;
192     int maxn    = 0;
193
194     for (int i = 0; i < iSizeB; i++)
195     {
196         n += piB[i];
197         sz += piB[i] * (piB[i] + 1) / 2;
198         upsz += piB[i] * piB[i];
199         maxn = std::max(maxn, piB[i]);
200     }
201
202     // optimal block size for dgels ????
203     int nb = 32;
204     int iSizeWork = (iSizeX + 2) * sz +
205                     upsz + 2 * n +
206                     std::max(std::max(iSizeX + sz * nb, 3 * maxn + maxn * (maxn + 1)), 3 * iSizeX);
207
208     double* pdblWork = new double[iSizeWork];
209     C2F(spf)(&iSizeX, &iSizeB, pDblF->get(), piB, pDblC->get(), pDblX->get(), pDblZ->get(),
210              pDblUl->get(), &dNu, &dAbstol, &dReltol, &dTv, &iIter, pdblWork, &iSizeWork, piWork, &iInfo);
211
212     delete pdblWork;
213     delete piWork;
214
215     if (iInfo < 0)
216     {
217         if (iInfo == -18)
218         {
219             Scierror(999, _("semi def fails.\n"));
220             return types::Function::Error;
221         }
222         else
223         {
224             Scierror(999, _("%s: Input argument %d of %s function has an illegal value.\n"), "semidef", -iInfo, "spf");
225             return types::Function::Error;
226         }
227     }
228
229     /*** return output arguments ***/
230     out.push_back(pDblX);
231
232     if (_iRetCount >= 2)
233     {
234         out.push_back(pDblZ);
235     }
236     else
237     {
238         delete pDblZ;
239     }
240
241     if (_iRetCount >= 3)
242     {
243         out.push_back(pDblUl);
244     }
245     else
246     {
247         delete pDblUl;
248     }
249
250     if (_iRetCount == 4)
251     {
252         types::Double* pDblOut = new types::Double(1, 2);
253         pDblOut->set(0, (double)iInfo);
254         pDblOut->set(1, (double)iIter);
255
256         out.push_back(pDblOut);
257     }
258
259     return types::Function::OK;
260 }
261