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 /*--------------------------------------------------------------------------*/
sci_semidef(types::typed_list & in,int _iRetCount,types::typed_list & out)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
131 // get C
132 if (in[4]->isDouble() == false)
133 {
134 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 5);
135 return types::Function::Error;
136 }
137
138 pDblC = in[4]->getAs<types::Double>();
139
140 if (pDblC->isComplex())
141 {
142 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 5);
143 return types::Function::Error;
144 }
145
146 if (pDblC->getRows()* pDblC->getCols() != pDblX->getRows()*pDblX->getCols())
147 {
148 Scierror(999, _("%s: Wrong size for input argument #%d: Input argument %d size expected.\n"), "semidef", 5, 1);
149 return types::Function::Error;
150 }
151
152 // get Option
153 if (in[5]->isDouble() == false)
154 {
155 Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "semidef", 6);
156 return types::Function::Error;
157 }
158
159 pDblOpt = in[5]->getAs<types::Double>();
160
161 if (pDblOpt->isComplex())
162 {
163 Scierror(999, _("%s: Wrong type for input argument #%d: A real vector expected.\n"), "semidef", 6);
164 return types::Function::Error;
165 }
166
167 if (pDblOpt->getRows() != 1 || pDblOpt->getCols() != 5)
168 {
169 Scierror(999, _("%s: Wrong size for input argument #%d: A vector of size %d x %d expected.\n"), "semidef", 6, 1, 5);
170 return types::Function::Error;
171 }
172
173 dNu = pDblOpt->get(0);
174 dAbstol = pDblOpt->get(1);
175 dReltol = pDblOpt->get(2);
176 dTv = pDblOpt->get(3);
177 iIter = (int)pDblOpt->get(4);
178
179 /*** perform operations ***/
180 types::Double* pDblUl = new types::Double(1, 2);
181 int* piWork = new int[iSizeX];
182
183 // compute size of working table
184 int n = 0;
185 int sz = 0;
186 int upsz = 0;
187 int maxn = 0;
188 piB = new int[iSizeB];
189 for (int i = 0; i < iSizeB; i++)
190 {
191 piB[i] = (int)pdblB[i];
192 }
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 delete[] piB;
215
216 if (iInfo < 0)
217 {
218 if (iInfo == -18)
219 {
220 Scierror(999, _("semi def fails.\n"));
221 return types::Function::Error;
222 }
223 else
224 {
225 Scierror(999, _("%s: Input argument %d of %s function has an illegal value.\n"), "semidef", -iInfo, "spf");
226 return types::Function::Error;
227 }
228 }
229
230 /*** return output arguments ***/
231 out.push_back(pDblX);
232
233 if (_iRetCount >= 2)
234 {
235 out.push_back(pDblZ);
236 }
237 else
238 {
239 delete pDblZ;
240 }
241
242 if (_iRetCount >= 3)
243 {
244 out.push_back(pDblUl);
245 }
246 else
247 {
248 delete pDblUl;
249 }
250
251 if (_iRetCount == 4)
252 {
253 types::Double* pDblOut = new types::Double(1, 2);
254 pDblOut->set(0, (double)iInfo);
255 pDblOut->set(1, (double)iIter);
256
257 out.push_back(pDblOut);
258 }
259
260 return types::Function::OK;
261 }
262
263