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