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);