10#define STAT_OUT plhs[0]
13static int DSDPPrintStats2(
DSDP,
void*);
14static int CountNonzeroMatrices(
double*,
int,
int,
int,
int*,
int*,
int*);
15static int CheckForConstantMat(
double*,
int,
int);
17static int printlevel=10;
18#define CHKERR(a) { if (a){mexWarnMsgTxt("DSDP Numerical Error"); } }
20#define mexErrMsgTxt(a); mexPrintf("Error: "); mexWarnMsgTxt(a); return;
35void mexFunction(
int nlhs, mxArray *plhs[],
36 int nrhs,
const mxArray *prhs[]){
38 mxArray *CA_cell_pr,*X_cell_pr;
39 const mxArray *OPTIONS_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;
46 int its,reuse=4,print_info=1,printsummary=0;
47 int ijnnz,spot,n,nn=0,nzmats,vecn;
49 int nsdpblocks=0,sdpblockj=0,sdpnmax=1,lpnmax=1,stat1=1,xmaker=0;
50 int sspot,nsubblocks,blockj;
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;
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"};
72 mexErrMsgTxt(
"Two input arguments required. See help for details. ");}
74 mexErrMsgTxt(
"Fewer input arguments required. See help for details. ");}
76 mexErrMsgTxt(
"Two output arguments required. See help for details. ");}
78 mexErrMsgTxt(
"Fewer output arguments required. See help for details. ");}
80 if (!mxIsDouble(B_IN) || mxIsSparse(B_IN)){
81 mexErrMsgTxt(
"DSDP: 1ST input must be a dense vector of doubles"); }
85 mexErrMsgTxt(
"DESP: 1ST input must be a column vector"); }
87 iscellCA = mxIsCell(CA_IN);
89 mexErrMsgTxt(
"DSDP: 2ND input must be a cell array"); }
93 mexErrMsgTxt(
"DSDP: dimension of 2ND cell array is p x 3");}
96 if(!mxIsStruct(PARS_IN)){
97 mexErrMsgTxt(
"3RD input `OPTIONS' should be a structure.");
102 if (!mxIsDouble(Y_IN) || mxIsSparse(Y_IN)){
103 mexErrMsgTxt(
"DSDP: 4TH input must be a dense vector of doubles"); }
106 if (m1 != nvars || n1 != nb){
107 mexErrMsgTxt(
"DSDP: dimensions of 1ST and 4TH input not compatible");}
110 mexErrMsgTxt(
"DSDP: Cannot read 4TH argument");}
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);
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);
127 subs[0] = j; subs[1] = 1;
128 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
129 CA_cell_pr = mxGetCell(CA_IN,index);
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);
141 subs[0] = j; subs[1] = 2;
142 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
143 CA_cell_pr = mxGetCell(CA_IN,index);
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."); }
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.");}
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;
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;}
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){
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);
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;
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 ){
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);
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.");}
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.");}
211 }
else if (strcmp(conetype,
"FIXED")==0){
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);
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);
234 mexPrintf(
"??? DSDP DATA ERROR: Cell: %d \n",j+1);
235 mexErrMsgTxt(
"DSDP: Third column in 2ND argument must have 1 row,");}
237 mexPrintf(
"??? DSDP DATA ERROR: Cell: %d \n",j+1);
238 mexErrMsgTxt(
"DSDP: Secord and third column must have same dimension,");}
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'. ");
248 if (X_OUT != NULL) mxDestroyArray(X_OUT) ;
249 X_OUT = mxCreateCellMatrix(mC,1);
251 if (Y_OUT != NULL) mxDestroyArray(Y_OUT) ;
252 Y_OUT = mxCreateDoubleMatrix(nvars, 1, mxREAL) ;
255 info = DSDPCreate(nvars,&dsdp); CHKERR(info);
256 info = DSDPCreateSDPCone(dsdp,nsdpblocks,&sdpcone); CHKERR(info);
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);}
264 for (j=0; j<mC; j++){
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){
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);
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;}
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");}
295 for (ii=0; ii<=nvars; ii++){
298 tnnz=0; spot=ajc[ii]; blockj=sdpblockj;
299 for (jj=0;jj<nsubblocks;jj++){
301 if (sdpnmax<n) sdpnmax=n;
304 info=SDPConeSetBlockSize(sdpcone,blockj,n); CHKERR(info);
305 info=SDPConeUsePackedFormat(sdpcone,blockj); CHKERR(info);
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;
312 for (tnnz2=tnnz+n*(n+1)/2,ijnnz=0;ijnnz<ajc[ii+1]-spot && air[spot+ijnnz]<tnnz2;ijnnz++){}
314 }
else if(ijnnz==n*(n+1)/2){
315 if (CheckForConstantMat(aval+spot,ijnnz,n)){
316 info=SDPConeSetConstantMat(sdpcone,blockj,i,n,aval[spot]); CHKERR(info);
318 info=SDPConeSetADenseVecMat(sdpcone,blockj,i,n,1.0,aval+spot,ijnnz); CHKERR(info);
321 info=SDPConeSetASparseVecMat(sdpcone,blockj,i,n,1.0,tnnz,air+spot,aval+spot,ijnnz); CHKERR(info);
323 tnnz+=n*(n+1)/2; spot+=ijnnz; blockj++;
326 sdpblockj=sdpblockj+nsubblocks;
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);
339 info=BConeSetUpperBound(bcone,i+1,aval[i]); CHKERR(info);
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);
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;
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);
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);
372 }
else if (strcmp(conetype,
"FIXED")==0){
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);
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");}
392 info=DSDPSetFixedVariables(dsdp,aval,aval2,xout,it1*it2); CHKERR(info);
393 for (i=0;i<it1*it2;i++){
405 for (i=0;i<nvars;i++){ info = DSDPSetY0(dsdp,i+1,y0[i]); CHKERR(info);}
408 reuse=(nvars-2)/sdpnmax;
409 if (nvars<50 && reuse==0) reuse=1;
410 if (reuse>=1) reuse++;
412 if (nvars<2000 && reuse>10) reuse=10;
413 if (reuse>12) reuse=12;
414 info=DSDPReuseMatrix(dsdp,reuse); CHKERR(info);
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);}
432 if(!mxIsStruct(PARS_IN)){
433 mexErrMsgTxt(
"Fifth Parameter `OPTIONS' should be a structure.");}
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;
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);}
523 info = DSDPSetMonitor(dsdp,DSDPPrintStats2,0); CHKERR(info);
525 info = DSDPSetup(dsdp);
527 mexErrMsgTxt(
"DSDP: Setup Error, Probably out of memory");}
530 info = DSDPSolve(dsdp); CHKERR(info);
531 info = DSDPStopReason(dsdp,&reason);
533 info=DSDPComputeX(dsdp);CHKERR(info);
535 info = DSDPGetSolutionType(dsdp,&pdfeasible); CHKERR(info);
537 if (printsummary){ DSDPEventLogSummary();}
540 mexErrMsgTxt(
"DSDP: Numerical error");}
543 mexErrMsgTxt(
"DSDP Terminated Due to Infeasible Starting Point\n");
544 }
else if (print_info){
547 mexPrintf(
"DSDP Converged. \n");
549 mexPrintf(
"DSDP Converged: Dual Objective exceeds its bound\n");
551 mexPrintf(
"DSDP Terminated Due to Small Steps\n");
553 mexPrintf(
"DSDP Terminated Due Maximum Number of Iterations\n");
555 mexPrintf(
"DSDP Terminated By User\n");
557 mexPrintf(
"DSDP Terminated Due to Infeasible Starting Point\n");
559 mexPrintf(
"DSDP Finished.\n");
563 mexPrintf(
"DSDP: Dual Unbounded, Primal Infeasible\n");
565 mexPrintf(
"DSDP: Primal Unbounded, Dual Infeasible\n");
570 info = DSDPGetY(dsdp,yout,nvars); CHKERR(info);
572 mexErrMsgTxt(
"DSDP: Numerical error");}
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);
591 STAT_FIELD = mxCreateString(
"Unbounded");
592 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
594 STAT_FIELD = mxCreateString(
"Infeasible");
595 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
597 STAT_FIELD = mxCreateString(
"PDFeasible");
598 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
601 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
602 stat=mxGetPr(STAT_FIELD); stat[0]=dobj;
603 mxSetField(STAT_OUT,0,fnames[1],STAT_FIELD);
605 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
606 stat=mxGetPr(STAT_FIELD); stat[0]=pobj;
607 mxSetField(STAT_OUT,0,fnames[2],STAT_FIELD);
609 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
610 stat=mxGetPr(STAT_FIELD); stat[0]=dobj;
611 mxSetField(STAT_OUT,0,fnames[3],STAT_FIELD);
614 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
615 stat=mxGetPr(STAT_FIELD); stat[0]=0;
616 mxSetField(STAT_OUT,0,fnames[4],STAT_FIELD);
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);
625 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
626 stat=mxGetPr(STAT_FIELD); stat[0]=0;
629 mxSetField(STAT_OUT,0,fnames[5],STAT_FIELD);
631 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
632 stat=mxGetPr(STAT_FIELD); stat[0]=-1;
636 mxSetField(STAT_OUT,0,fnames[5],STAT_FIELD);
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);
643 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
644 stat=mxGetPr(STAT_FIELD); stat[0]=dinf;
645 mxSetField(STAT_OUT,0,fnames[7],STAT_FIELD);
647 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
648 stat=mxGetPr(STAT_FIELD); stat[0]=mu;
649 mxSetField(STAT_OUT,0,fnames[8],STAT_FIELD);
651 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
652 stat=mxGetPr(STAT_FIELD); stat[0]=pstep;
653 mxSetField(STAT_OUT,0,fnames[9],STAT_FIELD);
655 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
656 stat=mxGetPr(STAT_FIELD); stat[0]=dstep;
657 mxSetField(STAT_OUT,0,fnames[10],STAT_FIELD);
659 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
660 stat=mxGetPr(STAT_FIELD); stat[0]=pnorm;
661 mxSetField(STAT_OUT,0,fnames[11],STAT_FIELD);
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);
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);
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);
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);
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);
690 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
691 stat=mxGetPr(STAT_FIELD);
693 mxSetField(STAT_OUT,0,fnames[17],STAT_FIELD);
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);
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);
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);
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);
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);
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);
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);
733 info = DSDPDestroy(dsdp); CHKERR(info);
741#define __FUNCT__ "CheckForConstantMat"
742static int CheckForConstantMat(
double*v,
int nnz,
int n){
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;}
753#define __FUNCT__ "CountNonzeroMatrices"
754 static int CountNonzeroMatrices(
double *blocksize,
int nblocks,
int block,
int nvars,
int *indd,
int*nnnz,
int *nnzmats){
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++){
762 while (indd[j]<marker1 && j< nnnz[i+1]){j++;}
763 if (j<nnnz[i+1] && indd[j]<marker2){ nzmats++;}
770 static int DSDPPrintStats2(
DSDP dsdp,
void* dummy){
773 double pobj,dobj,pstp=0,dstp,mu,res,pnorm,pinfeas;
776 if(printlevel<=0)
return(0);
778 info = DSDPStopReason(dsdp,&reason);
779 info = DSDPGetIts(dsdp,&iter);
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);
793 mexPrintf(
"Iter PP Objective DD Objective PInfeas DInfeas Nu StepLength Pnrm\n")
795 mexPrintf(
"----------------------------------------------------------------------------------------\n")
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);
801 mexPrintf(
" %1.0e \n",pnorm);
803 mexPrintf(
" %5.2f \n",pnorm);
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).
struct LPCone_C * LPCone
The LPCone object points to blocks of data that specify linear scalar inequality constraints.
DSDPTerminationReason
There are many reasons to terminate the solver.
DSDPSolutionType
Formulations (P) and (D) can be feasible and bounded, feasible and unbounded, or infeasible.
Internal structures for the DSDP solver.
Internal structure for semidefinite cone.