8 #define PARS_IN prhs[2] 10 #define STAT_OUT plhs[0] 13 static int DSDPPrintStats2(
DSDP,
void*);
14 static int CountNonzeroMatrices(
double*,
int,
int,
int,
int*,
int*,
int*);
15 static int CheckForConstantMat(
double*,
int,
int);
17 static int printlevel=10;
18 #define CHKERR(a) { if (a){mexWarnMsgTxt("DSDP Numerical Error"); } } 20 #define mexErrMsgTxt(a); mexPrintf("Error: "); mexWarnMsgTxt(a); return; 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) ;
256 info = DSDPCreateSDPCone(dsdp,nsdpblocks,&sdpcone); CHKERR(info);
260 if (!bval){ mexErrMsgTxt(
"DSDP: Problems with 1ST argument");}
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;
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);
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)){
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);
335 for (i=0;i<nvars;i++){
336 if (strcmp(conetype,
"LB")==0){
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");}
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");}
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;
417 info=
DSDPGetR(dsdp,&r0); 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);
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);
447 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"penalty") ){
448 penalty=mxGetScalar(OPTIONS_FIELD);
450 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"rho") ){
451 rho= mxGetScalar(OPTIONS_FIELD);
453 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"dynamicrho") ){
454 drho= (int)mxGetScalar(OPTIONS_FIELD);
456 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"zbar") ){
457 zbar= mxGetScalar(OPTIONS_FIELD);
459 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"dual_bound") ){
460 dbound= mxGetScalar(OPTIONS_FIELD);
462 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"reuse") ){
463 reuse= (int)mxGetScalar(OPTIONS_FIELD);
465 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"gaptol") ){
466 gaptol= mxGetScalar(OPTIONS_FIELD);
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);
479 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"inftol") ){
480 inftol= mxGetScalar(OPTIONS_FIELD);
482 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"infptol") ){
483 infptol= mxGetScalar(OPTIONS_FIELD);
485 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"pnormtol") ){
486 pnormtol= mxGetScalar(OPTIONS_FIELD);
488 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"boundy") ){
489 yhigh= fabs(mxGetScalar(OPTIONS_FIELD)); ylow=-yhigh;
491 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"r0") ){
492 r0= mxGetScalar(OPTIONS_FIELD);
494 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"mu0") ){
495 mu0=mxGetScalar(OPTIONS_FIELD);
497 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"maxtrustradius")){
498 maxtrust= mxGetScalar(OPTIONS_FIELD);
500 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"steptol") ){
501 steptol = mxGetScalar(OPTIONS_FIELD);
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") ){
513 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,
"xmaker") ){
514 xmaker = (int) mxGetScalar(OPTIONS_FIELD);}
523 info = DSDPSetMonitor(dsdp,DSDPPrintStats2,0); CHKERR(info);
527 mexErrMsgTxt(
"DSDP: Setup Error, Probably out of memory");}
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);
581 info=
DSDPGetR(dsdp,&dinf); 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);
667 mxSetField(STAT_OUT,0,fnames[12],STAT_FIELD);
669 STAT_FIELD = mxCreateDoubleMatrix(itmp, 1, mxREAL) ;
670 stat=mxGetPr(STAT_FIELD);
672 mxSetField(STAT_OUT,0,fnames[13],STAT_FIELD);
674 STAT_FIELD = mxCreateDoubleMatrix(6, 1, mxREAL) ;
675 stat=mxGetPr(STAT_FIELD);
677 mxSetField(STAT_OUT,0,fnames[14],STAT_FIELD);
679 STAT_FIELD = mxCreateDoubleMatrix(3, 1, mxREAL) ;
680 stat=mxGetPr(STAT_FIELD);
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);
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);
698 mxSetField(STAT_OUT,0,fnames[18],STAT_FIELD);
700 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
701 stat=mxGetPr(STAT_FIELD);
703 mxSetField(STAT_OUT,0,fnames[19],STAT_FIELD);
705 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
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);
713 mxSetField(STAT_OUT,0,fnames[21],STAT_FIELD);
716 STAT_FIELD = mxCreateDoubleMatrix(nvars+1, 1, mxREAL) ;
717 stat=mxGetPr(STAT_FIELD);
719 mxSetField(STAT_OUT,0,fnames[22],STAT_FIELD);
721 STAT_FIELD = mxCreateDoubleMatrix(nvars+1, 1, mxREAL) ;
722 stat=mxGetPr(STAT_FIELD);
724 mxSetField(STAT_OUT,0,fnames[23],STAT_FIELD);
726 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
727 stat=mxGetPr(STAT_FIELD);
729 mxSetField(STAT_OUT,0,fnames[24],STAT_FIELD);
741 #define __FUNCT__ "CheckForConstantMat" 742 static 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);
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);
int DSDPSetBarrierParameter(DSDP dsdp, double mu)
Set the current barrier parameter.
int DSDPCreate(int m, DSDP *dsdpnew)
Create a DSDP solver. FIRST DSDP routine!
int DSDPGetFinalErrors(DSDP dsdp, double err[6])
Copy six different error measurements into an array.
int DSDPGetDualBound(DSDP dsdp, double *dbound)
Get the termination parameter.
int DSDPDestroy(DSDP dsdp)
Free the internal data structures of the solver and the cones associated with it.
int DSDPGetStepTolerance(DSDP dsdp, double *steptol)
Get the current tolerance.
int DSDPGetStepLengths(DSDP dsdp, double *pstep, double *dstep)
Copy the step sizes in the current iteration.
int DSDPGetBarrierParameter(DSDP dsdp, double *mu)
Copy the current barrier parameter.
int DSDPGetDYMakeX(DSDP dsdp, double dy[], int m)
Copies the variables dy used to construct X into an array.
int DSDPGetPPObjective(DSDP dsdp, double *ppobj)
Copy the objective value (PP).
int DSDPAddObjectiveConstant(DSDP dsdp, double c)
Add a constant to the objective.
int DSDPGetRHistory(DSDP dsdp, double hist[], int length)
Copy a history of the infeasibility in (D) into an array.
int DSDPSetDualBound(DSDP dsdp, double dbound)
Terminate the solver if the objective value in (DD) is greater than this tolerance.
int DSDPGetPenaltyParameter(DSDP dsdp, double *Gamma)
Copy the penalty parameter Gamma.
int DSDPCreateBCone(DSDP dsdp, BCone *dspcone)
Create a new cone that represents bounds on the y variables.
int DSDPGetGapTolerance(DSDP dsdp, double *gaptol)
Get the termination tolerance.
int BConeSetUpperBound(BCone bcone, int vari, double ubound)
Set an upper bound on a variable y.
int DSDPSetZBar(DSDP dsdp, double ppobj)
Set an upper bound on the objective value at the solution.
int DSDPGetPTolerance(DSDP dsdp, double *inftol)
Copy the feasibility tolerance.
int DSDPSetDualObjective(DSDP dsdp, int i, double bi)
Set the objective vector b in (D).
int DSDPUseDynamicRho(DSDP dsdp, int yesorno)
Use a dynamic strategy to choose parameter rho.
int DSDPGetPnorm(DSDP dsdp, double *pnorm)
Copy the proximity of the solution to the central path.
int SDPConeSetXArray(SDPCone sdpcone, int blockj, int n, double xx[], int nn)
Provide an array for the SDPCone object can use to store dense matrices.
Internal structures for the DSDP solver.
int DSDPCreateLPCone(DSDP dsdp, LPCone *dspcone)
Create a new object for linear programs and scalar inequalities.
DSDPTerminationReason
There are many reasons to terminate the solver.
int DSDPSetGapTolerance(DSDP dsdp, double gaptol)
Terminate the solver when the relative duality gap is less than this tolerance.
int DSDPPrintData(DSDP dsdp, SDPCone sdpcone, LPCone lpcone)
Print data in SDPA format to a file named "output.sdpa".
int DSDPGetDataNorms(DSDP dsdp, double dnorm[3])
Copy the norms of the data C, A, and b into an array.
int BConeSetLowerBound(BCone bcone, int vari, double lbound)
Set a lower bound on a variable y.
The API to DSDP for those applications using DSDP as a subroutine library.
int DSDPSetRTolerance(DSDP dsdp, double inftol)
Classify (D) as feasible only if the variable r is less than this tolerance.
int DSDPSetPTolerance(DSDP dsdp, double inftol)
Classify (P) as feasible only if the infeasibility is less than this tolerance.
int SDPConeSetConstantMat(SDPCone sdpcone, int blockj, int vari, int n, double value)
Set a matrix whose elements are all the same.
int DSDPSetMaxTrustRadius(DSDP dsdp, double rad)
Set a maximum trust radius on the step direction.
int DSDPGetSolutionType(DSDP dsdp, DSDPSolutionType *pdfeasible)
Solutions can be bounded, infeasible, or unbounded.
int DSDPSetup(DSDP dsdp)
Set up data structures in the solver and the cones associated with it.
int DSDPReuseMatrix(DSDP dsdp, int rm)
Reuse the Hessian of the barrier function multiple times at each DSDP iteration.
int BConeAllocateBounds(BCone bcone, int nnz)
Set a surplus variable in constraint in (P).
int DSDPSetPenaltyParameter(DSDP dsdp, double Gamma)
Set the penalty parameter Gamma.
int DSDPGetReuseMatrix(DSDP dsdp, int *rm)
Copy this parameter.
int DSDPGetPInfeasibility(DSDP dsdp, double *pperror)
Copy the infeasibility in (P).
int DSDPGetMaxTrustRadius(DSDP dsdp, double *rad)
Copy the current radius of the trust region.
int DSDPGetDObjective(DSDP dsdp, double *dobj)
Copy the objective value (D).
int DSDPSetFixedVariables(DSDP dsdp, double vars[], double vals[], double xout[], int nvars)
Fix variable y to exact values.
int DSDPSetPotentialParameter(DSDP dsdp, double rho)
Set the potential parameter.
int DSDPStopReason(DSDP dsdp, DSDPTerminationReason *reason)
Copy the reason why the solver terminated.
int DSDPGetPObjective(DSDP dsdp, double *pobj)
Copy the objective value (P).
int DSDPSetY0(DSDP dsdp, int i, double yi0)
Set the initial values of variables y in (D).
int DSDPSetStepTolerance(DSDP dsdp, double steptol)
Terminate the solver if the step length in (DD) is below this tolerance.
int SDPConeUsePackedFormat(SDPCone sdpcone, int blockj)
Use packed symmetric format for the dense array.
int DSDPGetIts(DSDP dsdp, int *its)
Copy the current iteration number.
int DSDPSetPNormTolerance(DSDP dsdp, double ptol)
Terminate the solver when the relative duality gap is suffiently small and the PNorm is less than thi...
int DSDPGetMuMakeX(DSDP dsdp, double *mu)
Copies the value of mu used to construct X.
int DSDPGetYMaxNorm(DSDP dsdp, double *ynorm)
Copy the the infinity norm of the variables y.
int DSDPSolve(DSDP dsdp)
Apply DSDP to the problem.
int SDPConeSetBlockSize(SDPCone sdpcone, int blockj, int n)
Set the dimension of one block in the semidefinite cone.
Internal structure for semidefinite cone.
int DSDPGetY(DSDP dsdp, double y[], int m)
Copies the variables y into an array.
int SDPConeSetSparsity(SDPCone sdpcone, int blockj, int nnz)
Set the number of nonzero matrices in a block of the semidefinite cone.
int DSDPUsePenalty(DSDP dsdp, int yesorno)
Use penalty parameter to enforce feasibility.
int DSDPGetDDObjective(DSDP dsdp, double *ddobj)
Copy the objective value (DD).
int DSDPGetRTolerance(DSDP dsdp, double *inftol)
Copy the maximum infeasibility allowed (D).
int DSDPGetTraceX(DSDP dsdp, double *tracex)
Copy the trace of the variables X in (P).
int DSDPGetR(DSDP dsdp, double *res)
Copy the infeasibility in (D), or the variable r in (DD).
int DSDPComputeX(DSDP dsdp)
Compute the X variables.
int DSDPGetMaxIts(DSDP dsdp, int *its)
Copy the maximum number of iterations from the solver.
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
Call DSDP from the Matlab environment.
int SDPConeSetASparseVecMat(SDPCone sdpcone, int blockj, int vari, int n, double alpha, int ishift, const int ind[], const double val[], int nnz)
Set data matrix in a sparse format.
struct BCone_C * BCone
The BCone object points to lower and upper bounds on the variable y in (D).
int SDPConeSetADenseVecMat(SDPCone sdpcone, int blockj, int vari, int n, double alpha, double val[], int nnz)
Set a matrix in a dense format.
int DSDPSetYBounds(DSDP dsdp, double lbound, double ubound)
Bound the variables y.
int DSDPGetPotentialParameter(DSDP dsdp, double *rho)
Copy the potential parameter.
int DSDPGetYBounds(DSDP dsdp, double *lbound, double *ubound)
Copy the bounds on the variables y.
int LPConeSetData2(LPCone lpcone, int n, const int ik[], const int cols[], const double vals[])
Set data A and into the LP cone.
int DSDPSetR0(DSDP dsdp, double res)
Set an initial value for the variable r in (DD)
int DSDPSetMaxIts(DSDP dsdp, int its)
Terminate the solver after this number of iterations.
int DSDPGetGapHistory(DSDP dsdp, double hist[], int length)
Copy a history of the duality gap into an array.
int DSDPGetYMakeX(DSDP dsdp, double y[], int m)
Copies the variables y used to construct X into an array.
struct LPCone_C * LPCone
The LPCone object points to blocks of data that specify linear scalar inequality constraints.
DSDPSolutionType
Formulations (P) and (D) can be feasible and bounded, feasible and unbounded, or infeasible.