DSDP
dsdp.c
Go to the documentation of this file.
1#include "mex.h"
2#include "dsdp5.h"
3#include <math.h>
4
5/* #define K_IN prhs[0] */
6#define CA_IN prhs[1]
7#define B_IN prhs[0]
8#define PARS_IN prhs[2]
9#define Y_IN prhs[3]
10#define STAT_OUT plhs[0]
11#define Y_OUT plhs[1]
12#define X_OUT plhs[2]
13static int DSDPPrintStats2(DSDP, void*);
14static int CountNonzeroMatrices(double*,int,int,int,int*,int*,int*);
15static int CheckForConstantMat(double*,int,int);
16
17static int printlevel=10;
18#define CHKERR(a) { if (a){mexWarnMsgTxt("DSDP Numerical Error"); } }
19
20#define mexErrMsgTxt(a); mexPrintf("Error: "); mexWarnMsgTxt(a); return;
21
35void mexFunction(int nlhs, mxArray *plhs[],
36 int nrhs, const mxArray *prhs[]){
37
38 mxArray *CA_cell_pr,*X_cell_pr;
39 const mxArray *OPTIONS_FIELD;
40 mxArray *STAT_FIELD;
41 int i,j,k,ii,itmp,index,info;
42 int *air, *ajc, *air2, *ajc2, *str,*iptr;
43 int nvars,nb,mC,nC,m1,n1,m2,n2;
44 int nsubs=2, subs[2];
45 int iscellCA;
46 int its,reuse=4,print_info=1,printsummary=0;
47 int ijnnz,spot,n,nn=0,nzmats,vecn;
48 int it1,it2;
49 int nsdpblocks=0,sdpblockj=0,sdpnmax=1,lpnmax=1,stat1=1,xmaker=0;
50 int sspot,nsubblocks,blockj;
51 int jj,tnnz,tnnz2;
52 int maxit=1000,fastblas=1,rpos=0,drho=1,iloginfo=0,aggressive=0;
53 double penalty=1e8,rho=4,zbar=1e10,cc=0,r0=-1,mu0=-1,ylow,yhigh,gaptol=1e-6,pnormtol=1e30;
54 double maxtrust=1e30,steptol=0.01,inftol=1e-8,lpb=1.0,dbound=1e20,infptol=1e-4;
55 double dtmp,pstep,dstep,pnorm,mu;
56 double *blockn,datanorm[3];
57 double *aval,*aval2,*bval,*yout,*y0=0,*xout,*stat;
58 double pobj,dobj,dinf;
59 DSDP dsdp;
60 SDPCone sdpcone=0;
62 DSDPSolutionType pdfeasible;
63 LPCone lpcone=0;
64 BCone bcone=0;
65 char conetype[30];
66 int nfields=25;
67 const char *fnames[25]={"stype","obj","pobj","dobj","stopcode","termcode","iterates","r","mu",
68 "pstep","dstep","pnorm","gaphist","infeashist","errors",
69 "datanorm","ynorm","boundy","penalty","tracex","reuse","rho","xy","xdy","xmu"};
70
71 if (nrhs < 2){
72 mexErrMsgTxt("Two input arguments required. See help for details. ");}
73 if (nrhs > 4){
74 mexErrMsgTxt("Fewer input arguments required. See help for details. ");}
75 if (nlhs < 2){
76 mexErrMsgTxt("Two output arguments required. See help for details. ");}
77 if (nlhs > 3){
78 mexErrMsgTxt("Fewer output arguments required. See help for details. ");}
79
80 if (!mxIsDouble(B_IN) || mxIsSparse(B_IN)){
81 mexErrMsgTxt("DSDP: 1ST input must be a dense vector of doubles"); }
82 nvars = mxGetM(B_IN);
83 nb = mxGetN(B_IN);
84 if (nb > 1){
85 mexErrMsgTxt("DESP: 1ST input must be a column vector"); }
86
87 iscellCA = mxIsCell(CA_IN);
88 if (!iscellCA){
89 mexErrMsgTxt("DSDP: 2ND input must be a cell array"); }
90 mC = mxGetM(CA_IN);
91 nC = mxGetN(CA_IN);
92 if (nC != 3){
93 mexErrMsgTxt("DSDP: dimension of 2ND cell array is p x 3");}
94
95 if (nrhs >2){
96 if(!mxIsStruct(PARS_IN)){
97 mexErrMsgTxt("3RD input `OPTIONS' should be a structure.");
98 }
99 }
100
101 if (nrhs>3){
102 if (!mxIsDouble(Y_IN) || mxIsSparse(Y_IN)){
103 mexErrMsgTxt("DSDP: 4TH input must be a dense vector of doubles"); }
104 m1 = mxGetM(Y_IN);
105 n1 = mxGetN(Y_IN);
106 if (m1 != nvars || n1 != nb){
107 mexErrMsgTxt("DSDP: dimensions of 1ST and 4TH input not compatible");}
108 y0=mxGetPr(Y_IN);
109 if (!y0){
110 mexErrMsgTxt("DSDP: Cannot read 4TH argument");}
111 } else y0=0;
112
113
114 /* Check data */
115 for (j=0; j<mC; j++){
116 subs[0] = j; subs[1] = 0;
117 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
118 CA_cell_pr = mxGetCell(CA_IN,index);
119 if (!CA_cell_pr){
120 mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,1);
121 mexErrMsgTxt("DSDP: Empty Cell. Missing String"); }
122 if (!mxIsChar(CA_cell_pr)){
123 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
124 mexErrMsgTxt("DSDP: First column of cells in 2ND argument are a string to determine cone type."); }
125 mxGetString(CA_cell_pr,conetype,20);
126
127 subs[0] = j; subs[1] = 1;
128 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
129 CA_cell_pr = mxGetCell(CA_IN,index);
130 if (!CA_cell_pr){
131 mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,2);
132 mexErrMsgTxt("DSDP: Empty Cell. Provide dimension of block."); }
133 if (mxIsSparse(CA_cell_pr)){
134 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
135 mexErrMsgTxt("DSDP: Second column in 2ND argument must be a dense array of scalars that specify dimension."); }
136 if (!mxIsDouble(CA_cell_pr)){
137 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
138 mexErrMsgTxt("DSDP: Second column in 2ND argument must specify dimension."); }
139 aval=mxGetPr(CA_cell_pr);
140
141 subs[0] = j; subs[1] = 2;
142 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
143 CA_cell_pr = mxGetCell(CA_IN,index);
144 if (!CA_cell_pr){
145 mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,3);
146 mexErrMsgTxt("DSDP: Empty Cell. Provide sparse data matrix."); }
147 if (!mxIsSparse(CA_cell_pr) || !mxIsDouble(CA_cell_pr)){
148 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
149 mexErrMsgTxt("DSDP: Third column in 2ND argument must be a real sparse data matrix."); }
150
151 if (strcmp(conetype,"SDP")==0){
152 subs[0] = j; subs[1] = 1;
153 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
154 CA_cell_pr = mxGetCell(CA_IN,index);
155 aval=mxGetPr(CA_cell_pr);
156 it1 = mxGetM(CA_cell_pr);
157 it2 = mxGetN(CA_cell_pr);
158 if (it1!=1 && it2!=1){
159 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
160 mexErrMsgTxt("DSDP: Use a dense row vector in the second column in 2ND of the argument.");}
161
162 subs[0] = j; subs[1] = 2;
163 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
164 CA_cell_pr = mxGetCell(CA_IN,index);
165 m1 = mxGetN(CA_cell_pr)-1;
166 if (m1 != nvars){
167 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
168 mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have number of columns equal to number of variables+1.");}
169 vecn = mxGetM(CA_cell_pr);
170 for (tnnz=0,i=0;i<it1*it2;i++){n=(int)aval[i]; tnnz=tnnz+n*(n+1)/2;}
171 if ( tnnz != vecn){
172 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
173 mexErrMsgTxt("DSDP: Check Dimensions: The columns of A and C cannot be converted into square matrices");}
174 nsdpblocks=nsdpblocks+it1*it2;
175 } else if (strcmp(conetype,"LP")==0){
176
177 subs[0] = j; subs[1] = 1;
178 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
179 CA_cell_pr = mxGetCell(CA_IN,index);
180 it1 = mxGetM(CA_cell_pr);
181 it2 = mxGetN(CA_cell_pr);
182 /* THe sum of the dimensions should equal the number of constraints */
183
184 subs[0] = j; subs[1] = 2;
185 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
186 CA_cell_pr = mxGetCell(CA_IN,index);
187 if (!mxIsSparse(CA_cell_pr)){
188 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
189 mexErrMsgTxt("DSDP: Matrices in the third column in 2ND argument must be sparse."); }
190 m1 = mxGetN(CA_cell_pr)-1;
191 if (m1 != nvars){
192 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
193 mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have number of columns equal to number of variables+1.");}
194 } else if (strcmp(conetype,"LB")==0 || strcmp(conetype,"UB")==0 ){
195
196 subs[0] = j; subs[1] = 2;
197 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
198 CA_cell_pr = mxGetCell(CA_IN,index);
199 if (mxIsSparse(CA_cell_pr)){
200 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
201 mexErrMsgTxt("DSDP: Row vector in the third column in 2ND argument must be full."); }
202 it1 = mxGetM(CA_cell_pr);
203 it2 = mxGetN(CA_cell_pr);
204 if (it2 != 1){
205 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
206 mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have a single column of bounds.");}
207 if (it1 != nvars){
208 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
209 mexErrMsgTxt("DSDP: The column matrix in the third column in 2ND argument must contain a bound for each y variable.");}
210
211 } else if (strcmp(conetype,"FIXED")==0){
212 int dim1;
213 subs[0] = j; subs[1] = 1;
214 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
215 CA_cell_pr = mxGetCell(CA_IN,index);
216 if (mxIsSparse(CA_cell_pr)){
217 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
218 mexErrMsgTxt("DSDP: Vector in the third column in 2ND argument must be full."); }
219 it1 = mxGetM(CA_cell_pr);
220 it2 = mxGetN(CA_cell_pr);
221 dim1=it1*it2;
222 if (it1 != 1){
223 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
224 mexErrMsgTxt("DSDP: Third column in 2ND argument must have 1 row,");}
225 subs[0] = j; subs[1] = 2;
226 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
227 CA_cell_pr = mxGetCell(CA_IN,index);
228 if (mxIsSparse(CA_cell_pr)){
229 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
230 mexErrMsgTxt("DSDP: Vector in the third column in 2ND argument must be full."); }
231 it1 = mxGetM(CA_cell_pr);
232 it2 = mxGetN(CA_cell_pr);
233 if (it1 != 1){
234 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
235 mexErrMsgTxt("DSDP: Third column in 2ND argument must have 1 row,");}
236 if (it2 != dim1){
237 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
238 mexErrMsgTxt("DSDP: Secord and third column must have same dimension,");}
239 } else {
240 mexPrintf("??? DSDP DATA ERROR: Cell: %d, Conetype: %s \n",j+1,conetype);
241 mexErrMsgTxt("DSDP: Unknown Cone type in 2ND argument. Try 'SDP' or 'LP' or 'Bounds'. ");
242 }
243 }
244
245
246 /* Create output arrays */
247 if (nlhs>2){
248 if (X_OUT != NULL) mxDestroyArray(X_OUT) ;
249 X_OUT = mxCreateCellMatrix(mC,1);
250 }
251 if (Y_OUT != NULL) mxDestroyArray(Y_OUT) ;
252 Y_OUT = mxCreateDoubleMatrix(nvars, 1, mxREAL) ;
253
254 /* Create the Solver */
255 info = DSDPCreate(nvars,&dsdp); CHKERR(info);
256 info = DSDPCreateSDPCone(dsdp,nsdpblocks,&sdpcone); CHKERR(info);
257
258 /* Set Dual Objective Vector */
259 bval=mxGetPr(B_IN);
260 if (!bval){ mexErrMsgTxt("DSDP: Problems with 1ST argument");}
261 for (i=0;i<nvars;i++){info=DSDPSetDualObjective(dsdp,i+1,bval[i]); CHKERR(info);}
262
263 /* Set Matrix Data */
264 for (j=0; j<mC; j++){ /* Begin Block */
265 subs[0] = j; subs[1] = 0;
266 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
267 CA_cell_pr = mxGetCell(CA_IN,index);
268 mxGetString(CA_cell_pr,conetype,20);
269 if (strcmp(conetype,"SDP")==0){
270
271 subs[0] = j; subs[1] = 1;
272 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
273 CA_cell_pr = mxGetCell(CA_IN,index);
274 it1 = mxGetM(CA_cell_pr);
275 it2 = mxGetN(CA_cell_pr);
276 blockn=mxGetPr(CA_cell_pr);
277 nsubblocks=it1*it2;
278
279 subs[0] = j; subs[1] = 2;
280 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
281 CA_cell_pr = mxGetCell(CA_IN,index);
282 aval=mxGetPr(CA_cell_pr); air =mxGetIr(CA_cell_pr); ajc =mxGetJc(CA_cell_pr);
283 if (!aval||!air||!ajc)
284 { mexErrMsgTxt("DSDP: Problems with 2ND argument");}
285 for (tnnz=0,jj=0;jj<nsubblocks;jj++){n=(int)blockn[jj];tnnz+=n*(n+1)/2;}
286 if (nlhs>2){
287 subs[0] = j; subs[1] = 0;
288 index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
289 X_cell_pr = mxCreateDoubleMatrix(tnnz,1,mxREAL);
290 mxSetCell(X_OUT,index,X_cell_pr);
291 xout=mxGetPr(X_cell_pr);
292 if (tnnz>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
293 }
294
295 for (ii=0; ii<=nvars; ii++){ /* Begin Variable matrix constraints */
296 i=ii+1;
297 if (i==nvars+1) i=0;
298 tnnz=0; spot=ajc[ii]; blockj=sdpblockj;
299 for (jj=0;jj<nsubblocks;jj++){
300 n=(int)blockn[jj];
301 if (sdpnmax<n) sdpnmax=n;
302 if (ii==0){
303 nn+=n;
304 info=SDPConeSetBlockSize(sdpcone,blockj,n); CHKERR(info);
305 info=SDPConeUsePackedFormat(sdpcone,blockj); CHKERR(info);
306
307 if (nlhs>2){info=SDPConeSetXArray(sdpcone,blockj,n,xout+tnnz,(n*n+n)/2); CHKERR(info);}
308 info=CountNonzeroMatrices(blockn,nsubblocks,jj,nvars, air, ajc, &nzmats); CHKERR(info);
309 info=SDPConeSetSparsity(sdpcone,blockj,nzmats); CHKERR(info);
310 if (stat1<nzmats)stat1=nzmats;
311 }
312 for (tnnz2=tnnz+n*(n+1)/2,ijnnz=0;ijnnz<ajc[ii+1]-spot && air[spot+ijnnz]<tnnz2;ijnnz++){}
313 if ( ijnnz==0 ){ /* info=DSDPSetZeroMat(dsdp,sdpblockj,i,n); */
314 } else if(ijnnz==n*(n+1)/2){ /* check for dense matrix */
315 if (CheckForConstantMat(aval+spot,ijnnz,n)){
316 info=SDPConeSetConstantMat(sdpcone,blockj,i,n,aval[spot]); CHKERR(info);
317 } else {
318 info=SDPConeSetADenseVecMat(sdpcone,blockj,i,n,1.0,aval+spot,ijnnz); CHKERR(info);
319 }
320 } else {
321 info=SDPConeSetASparseVecMat(sdpcone,blockj,i,n,1.0,tnnz,air+spot,aval+spot,ijnnz); CHKERR(info);
322 }
323 tnnz+=n*(n+1)/2; spot+=ijnnz; blockj++;
324 }
325 } /* End Matrices in block */
326 sdpblockj=sdpblockj+nsubblocks;
327
328 } else if (strcmp(conetype,"LB")==0 || strcmp(conetype,"UB")==0 ){
329 subs[0] = j; subs[1] = 2;
330 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
331 CA_cell_pr = mxGetCell(CA_IN,index);
332 aval=mxGetPr(CA_cell_pr);
333 info=DSDPCreateBCone(dsdp,&bcone); CHKERR(info);
334 info=BConeAllocateBounds(bcone,nvars); CHKERR(info);
335 for (i=0;i<nvars;i++){
336 if (strcmp(conetype,"LB")==0){
337 info=BConeSetLowerBound(bcone,i+1,aval[i]); CHKERR(info);
338 } else {
339 info=BConeSetUpperBound(bcone,i+1,aval[i]); CHKERR(info);
340 }
341 }
342 if (nlhs>2){
343 subs[0] = j; subs[1] = 0;
344 index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
345 X_cell_pr = mxCreateDoubleMatrix(nvars,1,mxREAL);
346 mxSetCell(X_OUT,index,X_cell_pr);
347 aval2=mxGetPr(X_cell_pr);
348 info=BConeSetXArray(bcone,aval2,nvars); CHKERR(info);
349 }
350 } else if (strcmp(conetype,"LP")==0){
351 subs[0] = j; subs[1] = 2;
352 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
353 CA_cell_pr = mxGetCell(CA_IN,index);
354 n = mxGetM(CA_cell_pr);
355 if (lpnmax<n) lpnmax=n;
356 nn+=n;
357 aval=mxGetPr(CA_cell_pr); air =mxGetIr(CA_cell_pr); ajc =mxGetJc(CA_cell_pr);
358 if (!aval||!air||!ajc)
359 { mexErrMsgTxt("DSDP: Problems with 2ND argument");}
360 info=DSDPCreateLPCone(dsdp,&lpcone); CHKERR(info);
361 info=LPConeSetData2(lpcone,n,ajc,air,aval); CHKERR(info);
362 if (nlhs>2){
363 subs[0] = j; subs[1] = 0;
364 index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
365 X_cell_pr = mxCreateDoubleMatrix(n,1,mxREAL);
366 mxSetCell(X_OUT,index,X_cell_pr);
367 xout=mxGetPr(X_cell_pr);
368 if (n>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
369 info=LPConeSetXVec(lpcone,xout,n); CHKERR(info);
370 }
371
372 } else if (strcmp(conetype,"FIXED")==0){
373 int vari;
374 subs[0] = j; subs[1] = 1;
375 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
376 CA_cell_pr = mxGetCell(CA_IN,index);
377 aval=mxGetPr(CA_cell_pr);
378 it1 = mxGetM(CA_cell_pr);
379 it2 = mxGetN(CA_cell_pr);
380 subs[0] = j; subs[1] = 2;
381 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
382 CA_cell_pr = mxGetCell(CA_IN,index);
383 aval2=mxGetPr(CA_cell_pr);
384 if (nlhs>2){
385 subs[0] = j; subs[1] = 0;
386 index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
387 X_cell_pr = mxCreateDoubleMatrix(it1*it2,1,mxREAL);
388 mxSetCell(X_OUT,index,X_cell_pr);
389 xout=mxGetPr(X_cell_pr);
390 if (it1*it2>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
391 } else {xout=0;}
392 info=DSDPSetFixedVariables(dsdp,aval,aval2,xout,it1*it2); CHKERR(info);
393 for (i=0;i<it1*it2;i++){
394 /*
395 vari=(int)aval[i];
396 printf("FixedVariable %d to %4.4e\n",vari,aval2[i]);
397 info=DSDPSetFixedVariable(dsdp,vari,aval2[i]); CHKERR(info);
398 */
399 }
400 }
401 } /* End Block */
402
403 /* Set initial point */
404 if (y0){
405 for (i=0;i<nvars;i++){ info = DSDPSetY0(dsdp,i+1,y0[i]); CHKERR(info);}
406 }
407
408 reuse=(nvars-2)/sdpnmax;
409 if (nvars<50 && reuse==0) reuse=1;
410 if (reuse>=1) reuse++;
411 reuse=reuse*reuse;
412 if (nvars<2000 && reuse>10) reuse=10;
413 if (reuse>12) reuse=12;
414 info=DSDPReuseMatrix(dsdp,reuse); CHKERR(info);
415
416 info=DSDPGetPPObjective(dsdp,&zbar); CHKERR(info);
417 info=DSDPGetR(dsdp,&r0); CHKERR(info);
418 info=DSDPGetMaxIts(dsdp,&maxit); CHKERR(info);
419 info=DSDPGetPenaltyParameter(dsdp,&penalty); CHKERR(info);
420 info=DSDPGetPotentialParameter(dsdp,&rho); CHKERR(info);
421 info=DSDPGetDualBound(dsdp,&dbound); CHKERR(info);
422 info=DSDPGetGapTolerance(dsdp,&gaptol); CHKERR(info);
423 info=DSDPGetRTolerance(dsdp,&inftol); CHKERR(info);
424 info=DSDPGetBarrierParameter(dsdp,&mu0); CHKERR(info);
425 info=DSDPGetMaxTrustRadius(dsdp,&maxtrust); CHKERR(info);
426 info=DSDPGetStepTolerance(dsdp,&steptol); CHKERR(info);
427 info=DSDPGetPTolerance(dsdp,&infptol); CHKERR(info);
428 info=DSDPGetDataNorms(dsdp, datanorm); CHKERR(info);
429 info=DSDPGetYBounds(dsdp,&ylow,&yhigh); CHKERR(info);
430 if (datanorm[0]==0){info=DSDPSetYBounds(dsdp,-1.0,1.0);CHKERR(info);}
431 if (nrhs >2){
432 if(!mxIsStruct(PARS_IN)){
433 mexErrMsgTxt("Fifth Parameter `OPTIONS' should be a structure.");}
434
435 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"maxit") ){
436 maxit=(int) mxGetScalar(OPTIONS_FIELD);
437 info=DSDPSetMaxIts(dsdp,maxit); CHKERR(info);}
438 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"fastblas") ){
439 fastblas= (int) mxGetScalar(OPTIONS_FIELD);}
440 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"print") ){
441 print_info= (int) mxGetScalar(OPTIONS_FIELD);
442 printlevel=print_info;
443 }
444 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"bigM") ){
445 rpos=(int)mxGetScalar(OPTIONS_FIELD);
446 info=DSDPUsePenalty(dsdp,rpos); CHKERR(info);}
447 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"penalty") ){
448 penalty=mxGetScalar(OPTIONS_FIELD);
449 info=DSDPSetPenaltyParameter(dsdp,penalty); CHKERR(info);}
450 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"rho") ){
451 rho= mxGetScalar(OPTIONS_FIELD);
452 info=DSDPSetPotentialParameter(dsdp,rho); CHKERR(info);}
453 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dynamicrho") ){
454 drho= (int)mxGetScalar(OPTIONS_FIELD);
455 info=DSDPUseDynamicRho(dsdp,drho); CHKERR(info);}
456 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"zbar") ){
457 zbar= mxGetScalar(OPTIONS_FIELD);
458 info=DSDPSetZBar(dsdp,zbar); CHKERR(info);}
459 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dual_bound") ){
460 dbound= mxGetScalar(OPTIONS_FIELD);
461 info=DSDPSetDualBound(dsdp,dbound); CHKERR(info);}
462 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"reuse") ){
463 reuse= (int)mxGetScalar(OPTIONS_FIELD);
464 info=DSDPReuseMatrix(dsdp,reuse); CHKERR(info);}
465 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"gaptol") ){
466 gaptol= mxGetScalar(OPTIONS_FIELD);
467 info=DSDPSetGapTolerance(dsdp,gaptol);CHKERR(info);}
468 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"lp_barrier") ){
469 lpb= mxGetScalar(OPTIONS_FIELD);
470 if (lpb<0.1) lpb=0.1;
471 if (lpcone){info=LPConeScaleBarrier(lpcone,lpb); CHKERR(info);}}
472 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"lpb") ){
473 lpb= mxGetScalar(OPTIONS_FIELD);
474 if (lpb<0.1) lpb=0.1;
475 if (lpcone){info=LPConeScaleBarrier(lpcone,lpb); CHKERR(info);}}
476 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"cc") ){
477 cc= mxGetScalar(OPTIONS_FIELD);
478 info=DSDPAddObjectiveConstant(dsdp,cc); CHKERR(info);}
479 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"inftol") ){
480 inftol= mxGetScalar(OPTIONS_FIELD);
481 info=DSDPSetRTolerance(dsdp,inftol); CHKERR(info);}
482 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"infptol") ){
483 infptol= mxGetScalar(OPTIONS_FIELD);
484 info=DSDPSetPTolerance(dsdp,infptol);CHKERR(info);}
485 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"pnormtol") ){
486 pnormtol= mxGetScalar(OPTIONS_FIELD);
487 info=DSDPSetPNormTolerance(dsdp,pnormtol);CHKERR(info);}
488 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"boundy") ){
489 yhigh= fabs(mxGetScalar(OPTIONS_FIELD)); ylow=-yhigh;
490 info=DSDPSetYBounds(dsdp,ylow,yhigh);CHKERR(info);}
491 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"r0") ){
492 r0= mxGetScalar(OPTIONS_FIELD);
493 info=DSDPSetR0(dsdp,r0); CHKERR(info);}
494 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"mu0") ){
495 mu0=mxGetScalar(OPTIONS_FIELD);
496 info=DSDPSetBarrierParameter(dsdp,mu0);CHKERR(info);}
497 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"maxtrustradius")){
498 maxtrust= mxGetScalar(OPTIONS_FIELD);
499 info=DSDPSetMaxTrustRadius(dsdp,maxtrust); CHKERR(info);}
500 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"steptol") ){
501 steptol = mxGetScalar(OPTIONS_FIELD);
502 info=DSDPSetStepTolerance(dsdp,steptol); CHKERR(info);}
503 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dobjmin") ){
504 dtmp= mxGetScalar(OPTIONS_FIELD);
505 info = DSDPSetDualLowerBound(dsdp,dtmp);}
506 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dloginfo") ){
507 iloginfo= (int) mxGetScalar(OPTIONS_FIELD);
508 info=DSDPLogInfoAllow(iloginfo,0);}
509 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"logtime") ){
510 printsummary= (int) mxGetScalar(OPTIONS_FIELD);}
511 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"printproblem") ){
512 info=DSDPPrintData(dsdp,sdpcone,lpcone);}
513 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"xmaker") ){
514 xmaker = (int) mxGetScalar(OPTIONS_FIELD);}
515 /*
516 if( OPTIONS_FIELD = mxGetField(PARS_IN,0,"maxlanczos") ){
517 itmp= (int) mxGetScalar(OPTIONS_FIELD);
518 info=DSDPSetLanczosIterations(dsdp,itmp);}
519 */
520
521 }
522
523 info = DSDPSetMonitor(dsdp,DSDPPrintStats2,0); CHKERR(info);
524
525 info = DSDPSetup(dsdp);
526 if (info){
527 mexErrMsgTxt("DSDP: Setup Error, Probably out of memory");}
528
529
530 info = DSDPSolve(dsdp); CHKERR(info);
531 info = DSDPStopReason(dsdp,&reason);
532 if (reason!=DSDP_INFEASIBLE_START){
533 info=DSDPComputeX(dsdp);CHKERR(info);
534 }
535 info = DSDPGetSolutionType(dsdp,&pdfeasible); CHKERR(info);
536
537 if (printsummary){ DSDPEventLogSummary();}
538
539 if (info){
540 mexErrMsgTxt("DSDP: Numerical error");}
541
542 if ( reason == DSDP_INFEASIBLE_START){
543 mexErrMsgTxt("DSDP Terminated Due to Infeasible Starting Point\n");
544 } else if (print_info){
545
546 if (reason == DSDP_CONVERGED)
547 mexPrintf("DSDP Converged. \n");
548 else if ( reason == DSDP_UPPERBOUND )
549 mexPrintf("DSDP Converged: Dual Objective exceeds its bound\n");
550 else if ( reason == DSDP_SMALL_STEPS )
551 mexPrintf("DSDP Terminated Due to Small Steps\n");
552 else if ( reason == DSDP_MAX_IT)
553 mexPrintf("DSDP Terminated Due Maximum Number of Iterations\n");
554 else if ( reason == DSDP_USER_TERMINATION)
555 mexPrintf("DSDP Terminated By User\n");
556 else if ( reason == DSDP_INFEASIBLE_START)
557 mexPrintf("DSDP Terminated Due to Infeasible Starting Point\n");
558 else
559 mexPrintf("DSDP Finished.\n");
560 }
561
562 if (pdfeasible==DSDP_UNBOUNDED){
563 mexPrintf("DSDP: Dual Unbounded, Primal Infeasible\n");
564 } else if (pdfeasible==DSDP_INFEASIBLE){
565 mexPrintf("DSDP: Primal Unbounded, Dual Infeasible\n");
566 }
567
568 /* Set the dual solution */
569 yout=mxGetPr(Y_OUT);
570 info = DSDPGetY(dsdp,yout,nvars); CHKERR(info);
571 if (info){
572 mexErrMsgTxt("DSDP: Numerical error");}
573
574
575 /* Output statistics */
576 if (STAT_OUT != NULL) mxDestroyArray(STAT_OUT) ;
577 subs[0] = 1; subs[1] = 1;
578 STAT_OUT = mxCreateStructArray(2,subs,nfields,fnames);
579 info= DSDPGetDObjective(dsdp,&dobj); CHKERR(info);
580 info= DSDPGetPObjective(dsdp,&pobj); CHKERR(info);
581 info= DSDPGetR(dsdp,&dinf); CHKERR(info);
582 info= DSDPStopReason(dsdp,&reason); CHKERR(info);
583 info= DSDPGetIts(dsdp,&its); CHKERR(info);
584 info= DSDPGetBarrierParameter(dsdp,&mu); CHKERR(info);
585 info= DSDPGetStepLengths(dsdp,&pstep,&dstep); CHKERR(info);
586 info= DSDPGetPnorm(dsdp,&pnorm); CHKERR(info);
587 info= DSDPGetBarrierParameter(dsdp,&mu); CHKERR(info);
588 info= DSDPGetYBounds(dsdp,&ylow,&yhigh); CHKERR(info);
589
590 if (pdfeasible==DSDP_UNBOUNDED){
591 STAT_FIELD = mxCreateString("Unbounded");
592 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
593 } else if (pdfeasible==DSDP_INFEASIBLE){
594 STAT_FIELD = mxCreateString("Infeasible");
595 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
596 } else {
597 STAT_FIELD = mxCreateString("PDFeasible");
598 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
599 }
600
601 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
602 stat=mxGetPr(STAT_FIELD); stat[0]=dobj;
603 mxSetField(STAT_OUT,0,fnames[1],STAT_FIELD);
604
605 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
606 stat=mxGetPr(STAT_FIELD); stat[0]=pobj;
607 mxSetField(STAT_OUT,0,fnames[2],STAT_FIELD);
608
609 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
610 stat=mxGetPr(STAT_FIELD); stat[0]=dobj;
611 mxSetField(STAT_OUT,0,fnames[3],STAT_FIELD);
612
613 if (reason==DSDP_CONVERGED){
614 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
615 stat=mxGetPr(STAT_FIELD); stat[0]=0;
616 mxSetField(STAT_OUT,0,fnames[4],STAT_FIELD);
617 } else {
618 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
619 stat=mxGetPr(STAT_FIELD); stat[0]=(double)reason;
620 if (stat[0]==0) stat[0]=-1;
621 mxSetField(STAT_OUT,0,fnames[4],STAT_FIELD);
622 }
623
624 if (reason==DSDP_CONVERGED){ /* For YALMIP */
625 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
626 stat=mxGetPr(STAT_FIELD); stat[0]=0;
627 if (pdfeasible==DSDP_UNBOUNDED){ stat[0]=1;}
628 if (pdfeasible==DSDP_INFEASIBLE){ stat[0]=2;}
629 mxSetField(STAT_OUT,0,fnames[5],STAT_FIELD);
630 } else {
631 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
632 stat=mxGetPr(STAT_FIELD); stat[0]=-1;
633 if (reason==DSDP_INFEASIBLE_START)stat[0]=-6;
634 if (reason==DSDP_UPPERBOUND)stat[0]=-5;
635 if (reason==DSDP_MAX_IT)stat[0]=-3;
636 mxSetField(STAT_OUT,0,fnames[5],STAT_FIELD);
637 }
638
639 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
640 stat=mxGetPr(STAT_FIELD); stat[0]=(double) its;
641 mxSetField(STAT_OUT,0,fnames[6],STAT_FIELD);
642
643 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
644 stat=mxGetPr(STAT_FIELD); stat[0]=dinf;
645 mxSetField(STAT_OUT,0,fnames[7],STAT_FIELD);
646
647 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
648 stat=mxGetPr(STAT_FIELD); stat[0]=mu;
649 mxSetField(STAT_OUT,0,fnames[8],STAT_FIELD);
650
651 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
652 stat=mxGetPr(STAT_FIELD); stat[0]=pstep;
653 mxSetField(STAT_OUT,0,fnames[9],STAT_FIELD);
654
655 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
656 stat=mxGetPr(STAT_FIELD); stat[0]=dstep;
657 mxSetField(STAT_OUT,0,fnames[10],STAT_FIELD);
658
659 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
660 stat=mxGetPr(STAT_FIELD); stat[0]=pnorm;
661 mxSetField(STAT_OUT,0,fnames[11],STAT_FIELD);
662
663 itmp=100; if (its < itmp) itmp=its;
664 STAT_FIELD = mxCreateDoubleMatrix(itmp, 1, mxREAL) ;
665 stat=mxGetPr(STAT_FIELD);
666 info= DSDPGetGapHistory(dsdp,stat,itmp); CHKERR(info);
667 mxSetField(STAT_OUT,0,fnames[12],STAT_FIELD);
668
669 STAT_FIELD = mxCreateDoubleMatrix(itmp, 1, mxREAL) ;
670 stat=mxGetPr(STAT_FIELD);
671 info= DSDPGetRHistory(dsdp,stat,itmp); CHKERR(info);
672 mxSetField(STAT_OUT,0,fnames[13],STAT_FIELD);
673
674 STAT_FIELD = mxCreateDoubleMatrix(6, 1, mxREAL) ;
675 stat=mxGetPr(STAT_FIELD);
676 info = DSDPGetFinalErrors(dsdp,stat); CHKERR(info);
677 mxSetField(STAT_OUT,0,fnames[14],STAT_FIELD);
678
679 STAT_FIELD = mxCreateDoubleMatrix(3, 1, mxREAL) ;
680 stat=mxGetPr(STAT_FIELD);
681 info = DSDPGetDataNorms(dsdp,stat); CHKERR(info);
682 dtmp=stat[0];stat[0]=stat[1];stat[1]=dtmp;
683 mxSetField(STAT_OUT,0,fnames[15],STAT_FIELD);
684
685 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
686 stat=mxGetPr(STAT_FIELD);
687 info = DSDPGetYMaxNorm(dsdp,stat); CHKERR(info);
688 mxSetField(STAT_OUT,0,fnames[16],STAT_FIELD);
689
690 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
691 stat=mxGetPr(STAT_FIELD);
692 stat[0]=yhigh;
693 mxSetField(STAT_OUT,0,fnames[17],STAT_FIELD);
694
695 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
696 stat=mxGetPr(STAT_FIELD);
697 info = DSDPGetPenaltyParameter(dsdp,stat); CHKERR(info);
698 mxSetField(STAT_OUT,0,fnames[18],STAT_FIELD);
699
700 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
701 stat=mxGetPr(STAT_FIELD);
702 info = DSDPGetTraceX(dsdp,stat); CHKERR(info);
703 mxSetField(STAT_OUT,0,fnames[19],STAT_FIELD);
704
705 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
706 info = DSDPGetReuseMatrix(dsdp,&its); CHKERR(info);
707 stat=mxGetPr(STAT_FIELD); stat[0]=(double) its;
708 mxSetField(STAT_OUT,0,fnames[20],STAT_FIELD);
709
710 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
711 stat=mxGetPr(STAT_FIELD);
712 info = DSDPGetPotentialParameter(dsdp,stat); CHKERR(info);
713 mxSetField(STAT_OUT,0,fnames[21],STAT_FIELD);
714
715 if (xmaker){
716 STAT_FIELD = mxCreateDoubleMatrix(nvars+1, 1, mxREAL) ;
717 stat=mxGetPr(STAT_FIELD);
718 info = DSDPGetYMakeX(dsdp,stat,nvars+1); CHKERR(info);
719 mxSetField(STAT_OUT,0,fnames[22],STAT_FIELD);
720
721 STAT_FIELD = mxCreateDoubleMatrix(nvars+1, 1, mxREAL) ;
722 stat=mxGetPr(STAT_FIELD);
723 info = DSDPGetDYMakeX(dsdp,stat,nvars+1); CHKERR(info);
724 mxSetField(STAT_OUT,0,fnames[23],STAT_FIELD);
725
726 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
727 stat=mxGetPr(STAT_FIELD);
728 info = DSDPGetMuMakeX(dsdp,stat); CHKERR(info);
729 mxSetField(STAT_OUT,0,fnames[24],STAT_FIELD);
730 }
731 /* Free internal data structure */
732
733 info = DSDPDestroy(dsdp); CHKERR(info);
734
735 return;
736} /* main */
737
738
739
740#undef __FUNCT__
741#define __FUNCT__ "CheckForConstantMat"
742static int CheckForConstantMat(double*v,int nnz, int n){
743 int i;double vv;
744 if (n<=1){ return 0; }
745 if (nnz!=(n*n+n)/2){ return 0; }
746 for (vv=v[0],i=1;i<nnz;i++){
747 if (v[i]!=vv){ return 0;}
748 }
749 return 1;
750}
751
752#undef __FUNCT__
753#define __FUNCT__ "CountNonzeroMatrices"
754 static int CountNonzeroMatrices(double *blocksize, int nblocks, int block, int nvars, int *indd, int*nnnz, int *nnzmats){
755 int i,j,n,nzmats=0;
756 int marker1=0,marker2=0;
757 for (i=0;i<block;i++){n=(int)blocksize[i]; marker1=marker1+n*(n+1)/2;}
758 n=(int)blocksize[block];
759 marker2=marker1+n*(n+1)/2;
760 for (i=0;i<nvars;i++){
761 j=nnnz[i];
762 while (indd[j]<marker1 && j< nnnz[i+1]){j++;}
763 if (j<nnnz[i+1] && indd[j]<marker2){ nzmats++;}
764 }
765 *nnzmats=nzmats;
766 return 0;
767}
768
769/* ---------------------------------------------------------- */
770 static int DSDPPrintStats2(DSDP dsdp, void* dummy){
771
772 int iter,info;
773 double pobj,dobj,pstp=0,dstp,mu,res,pnorm,pinfeas;
775
776 if(printlevel<=0) return(0);
777
778 info = DSDPStopReason(dsdp,&reason);
779 info = DSDPGetIts(dsdp,&iter);
780
781 if( (reason!=CONTINUE_ITERATING) || ((iter % printlevel)==0)){
782
783 info = DSDPGetDDObjective(dsdp,&dobj);
784 info = DSDPGetPPObjective(dsdp,&pobj);
785 info = DSDPGetR(dsdp,&res);
786 info = DSDPGetPInfeasibility(dsdp,&pinfeas);
787 info = DSDPGetStepLengths(dsdp,&pstp,&dstp);
788 info = DSDPGetBarrierParameter(dsdp,&mu);
789 info = DSDPGetPnorm(dsdp,&pnorm);
790
791
792 if (iter==0){
793 mexPrintf("Iter PP Objective DD Objective PInfeas DInfeas Nu StepLength Pnrm\n")
794 ;
795 mexPrintf("----------------------------------------------------------------------------------------\n")
796 ;
797 }
798 mexPrintf("%-3d %16.8e %16.8e %9.1e %9.1e %9.1e",iter,pobj,dobj,pinfeas,res,mu);
799 mexPrintf(" %4.2f %4.2f",pstp,dstp);
800 if (pnorm>1.0e3){
801 mexPrintf(" %1.0e \n",pnorm);
802 } else {
803 mexPrintf(" %5.2f \n",pnorm);
804 }
805 fflush(NULL);
806 }
807 return(0);
808}
The API to DSDP for those applications using DSDP as a subroutine library.
struct BCone_C * BCone
The BCone object points to lower and upper bounds on the variable y in (D).
Definition dsdp5.h:28
struct LPCone_C * LPCone
The LPCone object points to blocks of data that specify linear scalar inequality constraints.
Definition dsdp5.h:27
DSDPTerminationReason
There are many reasons to terminate the solver.
@ DSDP_USER_TERMINATION
@ CONTINUE_ITERATING
@ DSDP_MAX_IT
@ DSDP_UPPERBOUND
@ DSDP_INFEASIBLE_START
@ DSDP_CONVERGED
@ DSDP_SMALL_STEPS
DSDPSolutionType
Formulations (P) and (D) can be feasible and bounded, feasible and unbounded, or infeasible.
@ DSDP_UNBOUNDED
@ DSDP_INFEASIBLE
Internal structures for the DSDP solver.
Definition dsdp.h:65
Internal structure for semidefinite cone.
Definition dsdpsdp.h:80