00001 #include "mex.h"
00002 #include "dsdp5.h"
00003 #include <math.h>
00004
00005
00006 #define CA_IN prhs[1]
00007 #define B_IN prhs[0]
00008 #define PARS_IN prhs[2]
00009 #define Y_IN prhs[3]
00010 #define STAT_OUT plhs[0]
00011 #define Y_OUT plhs[1]
00012 #define X_OUT plhs[2]
00013 static int DSDPPrintStats2(DSDP, void*);
00014 static int CountNonzeroMatrices(double*,int,int,int,int*,int*,int*);
00015 static int CheckForConstantMat(double*,int,int);
00016
00017 static int printlevel=10;
00018 #define CHKERR(a) { if (a){mexWarnMsgTxt("DSDP Numerical Error"); } }
00019
00020 #define mexErrMsgTxt(a); mexPrintf("Error: "); mexWarnMsgTxt(a); return;
00021
00035 void mexFunction(int nlhs, mxArray *plhs[],
00036 int nrhs, const mxArray *prhs[]){
00037
00038 mxArray *CA_cell_pr,*X_cell_pr;
00039 const mxArray *OPTIONS_FIELD;
00040 mxArray *STAT_FIELD;
00041 int i,j,k,ii,itmp,index,info;
00042 int *air, *ajc, *air2, *ajc2, *str,*iptr;
00043 int nvars,nb,mC,nC,m1,n1,m2,n2;
00044 int nsubs=2, subs[2];
00045 int iscellCA;
00046 int its,reuse=4,print_info=1,printsummary=0;
00047 int ijnnz,spot,n,nn=0,nzmats,vecn;
00048 int it1,it2;
00049 int nsdpblocks=0,sdpblockj=0,sdpnmax=1,lpnmax=1,stat1=1,xmaker=0;
00050 int sspot,nsubblocks,blockj;
00051 int jj,tnnz,tnnz2;
00052 int maxit=1000,fastblas=1,rpos=0,drho=1,iloginfo=0,aggressive=0;
00053 double penalty=1e8,rho=4,zbar=1e10,cc=0,r0=-1,mu0=-1,ylow,yhigh,gaptol=1e-6,pnormtol=1e30;
00054 double maxtrust=1e30,steptol=0.01,inftol=1e-8,lpb=1.0,dbound=1e20,infptol=1e-4;
00055 double dtmp,pstep,dstep,pnorm,mu;
00056 double *blockn,datanorm[3];
00057 double *aval,*aval2,*bval,*yout,*y0=0,*xout,*stat;
00058 double pobj,dobj,dinf;
00059 DSDP dsdp;
00060 SDPCone sdpcone=0;
00061 DSDPTerminationReason reason;
00062 DSDPSolutionType pdfeasible;
00063 LPCone lpcone=0;
00064 BCone bcone=0;
00065 char conetype[30];
00066 int nfields=25;
00067 const char *fnames[25]={"stype","obj","pobj","dobj","stopcode","termcode","iterates","r","mu",
00068 "pstep","dstep","pnorm","gaphist","infeashist","errors",
00069 "datanorm","ynorm","boundy","penalty","tracex","reuse","rho","xy","xdy","xmu"};
00070
00071 if (nrhs < 2){
00072 mexErrMsgTxt("Two input arguments required. See help for details. ");}
00073 if (nrhs > 4){
00074 mexErrMsgTxt("Fewer input arguments required. See help for details. ");}
00075 if (nlhs < 2){
00076 mexErrMsgTxt("Two output arguments required. See help for details. ");}
00077 if (nlhs > 3){
00078 mexErrMsgTxt("Fewer output arguments required. See help for details. ");}
00079
00080 if (!mxIsDouble(B_IN) || mxIsSparse(B_IN)){
00081 mexErrMsgTxt("DSDP: 1ST input must be a dense vector of doubles"); }
00082 nvars = mxGetM(B_IN);
00083 nb = mxGetN(B_IN);
00084 if (nb > 1){
00085 mexErrMsgTxt("DESP: 1ST input must be a column vector"); }
00086
00087 iscellCA = mxIsCell(CA_IN);
00088 if (!iscellCA){
00089 mexErrMsgTxt("DSDP: 2ND input must be a cell array"); }
00090 mC = mxGetM(CA_IN);
00091 nC = mxGetN(CA_IN);
00092 if (nC != 3){
00093 mexErrMsgTxt("DSDP: dimension of 2ND cell array is p x 3");}
00094
00095 if (nrhs >2){
00096 if(!mxIsStruct(PARS_IN)){
00097 mexErrMsgTxt("3RD input `OPTIONS' should be a structure.");
00098 }
00099 }
00100
00101 if (nrhs>3){
00102 if (!mxIsDouble(Y_IN) || mxIsSparse(Y_IN)){
00103 mexErrMsgTxt("DSDP: 4TH input must be a dense vector of doubles"); }
00104 m1 = mxGetM(Y_IN);
00105 n1 = mxGetN(Y_IN);
00106 if (m1 != nvars || n1 != nb){
00107 mexErrMsgTxt("DSDP: dimensions of 1ST and 4TH input not compatible");}
00108 y0=mxGetPr(Y_IN);
00109 if (!y0){
00110 mexErrMsgTxt("DSDP: Cannot read 4TH argument");}
00111 } else y0=0;
00112
00113
00114
00115 for (j=0; j<mC; j++){
00116 subs[0] = j; subs[1] = 0;
00117 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00118 CA_cell_pr = mxGetCell(CA_IN,index);
00119 if (!CA_cell_pr){
00120 mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,1);
00121 mexErrMsgTxt("DSDP: Empty Cell. Missing String"); }
00122 if (!mxIsChar(CA_cell_pr)){
00123 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00124 mexErrMsgTxt("DSDP: First column of cells in 2ND argument are a string to determine cone type."); }
00125 mxGetString(CA_cell_pr,conetype,20);
00126
00127 subs[0] = j; subs[1] = 1;
00128 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00129 CA_cell_pr = mxGetCell(CA_IN,index);
00130 if (!CA_cell_pr){
00131 mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,2);
00132 mexErrMsgTxt("DSDP: Empty Cell. Provide dimension of block."); }
00133 if (mxIsSparse(CA_cell_pr)){
00134 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00135 mexErrMsgTxt("DSDP: Second column in 2ND argument must be a dense array of scalars that specify dimension."); }
00136 if (!mxIsDouble(CA_cell_pr)){
00137 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00138 mexErrMsgTxt("DSDP: Second column in 2ND argument must specify dimension."); }
00139 aval=mxGetPr(CA_cell_pr);
00140
00141 subs[0] = j; subs[1] = 2;
00142 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00143 CA_cell_pr = mxGetCell(CA_IN,index);
00144 if (!CA_cell_pr){
00145 mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,3);
00146 mexErrMsgTxt("DSDP: Empty Cell. Provide sparse data matrix."); }
00147 if (!mxIsSparse(CA_cell_pr) || !mxIsDouble(CA_cell_pr)){
00148 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00149 mexErrMsgTxt("DSDP: Third column in 2ND argument must be a real sparse data matrix."); }
00150
00151 if (strcmp(conetype,"SDP")==0){
00152 subs[0] = j; subs[1] = 1;
00153 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00154 CA_cell_pr = mxGetCell(CA_IN,index);
00155 aval=mxGetPr(CA_cell_pr);
00156 it1 = mxGetM(CA_cell_pr);
00157 it2 = mxGetN(CA_cell_pr);
00158 if (it1!=1 && it2!=1){
00159 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00160 mexErrMsgTxt("DSDP: Use a dense row vector in the second column in 2ND of the argument.");}
00161
00162 subs[0] = j; subs[1] = 2;
00163 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00164 CA_cell_pr = mxGetCell(CA_IN,index);
00165 m1 = mxGetN(CA_cell_pr)-1;
00166 if (m1 != nvars){
00167 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00168 mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have number of columns equal to number of variables+1.");}
00169 vecn = mxGetM(CA_cell_pr);
00170 for (tnnz=0,i=0;i<it1*it2;i++){n=(int)aval[i]; tnnz=tnnz+n*(n+1)/2;}
00171 if ( tnnz != vecn){
00172 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00173 mexErrMsgTxt("DSDP: Check Dimensions: The columns of A and C cannot be converted into square matrices");}
00174 nsdpblocks=nsdpblocks+it1*it2;
00175 } else if (strcmp(conetype,"LP")==0){
00176
00177 subs[0] = j; subs[1] = 1;
00178 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00179 CA_cell_pr = mxGetCell(CA_IN,index);
00180 it1 = mxGetM(CA_cell_pr);
00181 it2 = mxGetN(CA_cell_pr);
00182
00183
00184 subs[0] = j; subs[1] = 2;
00185 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00186 CA_cell_pr = mxGetCell(CA_IN,index);
00187 if (!mxIsSparse(CA_cell_pr)){
00188 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00189 mexErrMsgTxt("DSDP: Matrices in the third column in 2ND argument must be sparse."); }
00190 m1 = mxGetN(CA_cell_pr)-1;
00191 if (m1 != nvars){
00192 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00193 mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have number of columns equal to number of variables+1.");}
00194 } else if (strcmp(conetype,"LB")==0 || strcmp(conetype,"UB")==0 ){
00195
00196 subs[0] = j; subs[1] = 2;
00197 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00198 CA_cell_pr = mxGetCell(CA_IN,index);
00199 if (mxIsSparse(CA_cell_pr)){
00200 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00201 mexErrMsgTxt("DSDP: Row vector in the third column in 2ND argument must be full."); }
00202 it1 = mxGetM(CA_cell_pr);
00203 it2 = mxGetN(CA_cell_pr);
00204 if (it2 != 1){
00205 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00206 mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have a single column of bounds.");}
00207 if (it1 != nvars){
00208 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00209 mexErrMsgTxt("DSDP: The column matrix in the third column in 2ND argument must contain a bound for each y variable.");}
00210
00211 } else if (strcmp(conetype,"FIXED")==0){
00212 int dim1;
00213 subs[0] = j; subs[1] = 1;
00214 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00215 CA_cell_pr = mxGetCell(CA_IN,index);
00216 if (mxIsSparse(CA_cell_pr)){
00217 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00218 mexErrMsgTxt("DSDP: Vector in the third column in 2ND argument must be full."); }
00219 it1 = mxGetM(CA_cell_pr);
00220 it2 = mxGetN(CA_cell_pr);
00221 dim1=it1*it2;
00222 if (it1 != 1){
00223 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00224 mexErrMsgTxt("DSDP: Third column in 2ND argument must have 1 row,");}
00225 subs[0] = j; subs[1] = 2;
00226 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00227 CA_cell_pr = mxGetCell(CA_IN,index);
00228 if (mxIsSparse(CA_cell_pr)){
00229 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00230 mexErrMsgTxt("DSDP: Vector in the third column in 2ND argument must be full."); }
00231 it1 = mxGetM(CA_cell_pr);
00232 it2 = mxGetN(CA_cell_pr);
00233 if (it1 != 1){
00234 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00235 mexErrMsgTxt("DSDP: Third column in 2ND argument must have 1 row,");}
00236 if (it2 != dim1){
00237 mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00238 mexErrMsgTxt("DSDP: Secord and third column must have same dimension,");}
00239 } else {
00240 mexPrintf("??? DSDP DATA ERROR: Cell: %d, Conetype: %s \n",j+1,conetype);
00241 mexErrMsgTxt("DSDP: Unknown Cone type in 2ND argument. Try 'SDP' or 'LP' or 'Bounds'. ");
00242 }
00243 }
00244
00245
00246
00247 if (nlhs>2){
00248 if (X_OUT != NULL) mxDestroyArray(X_OUT) ;
00249 X_OUT = mxCreateCellMatrix(mC,1);
00250 }
00251 if (Y_OUT != NULL) mxDestroyArray(Y_OUT) ;
00252 Y_OUT = mxCreateDoubleMatrix(nvars, 1, mxREAL) ;
00253
00254
00255 info = DSDPCreate(nvars,&dsdp); CHKERR(info);
00256 info = DSDPCreateSDPCone(dsdp,nsdpblocks,&sdpcone); CHKERR(info);
00257
00258
00259 bval=mxGetPr(B_IN);
00260 if (!bval){ mexErrMsgTxt("DSDP: Problems with 1ST argument");}
00261 for (i=0;i<nvars;i++){info=DSDPSetDualObjective(dsdp,i+1,bval[i]); CHKERR(info);}
00262
00263
00264 for (j=0; j<mC; j++){
00265 subs[0] = j; subs[1] = 0;
00266 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00267 CA_cell_pr = mxGetCell(CA_IN,index);
00268 mxGetString(CA_cell_pr,conetype,20);
00269 if (strcmp(conetype,"SDP")==0){
00270
00271 subs[0] = j; subs[1] = 1;
00272 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00273 CA_cell_pr = mxGetCell(CA_IN,index);
00274 it1 = mxGetM(CA_cell_pr);
00275 it2 = mxGetN(CA_cell_pr);
00276 blockn=mxGetPr(CA_cell_pr);
00277 nsubblocks=it1*it2;
00278
00279 subs[0] = j; subs[1] = 2;
00280 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00281 CA_cell_pr = mxGetCell(CA_IN,index);
00282 aval=mxGetPr(CA_cell_pr); air =mxGetIr(CA_cell_pr); ajc =mxGetJc(CA_cell_pr);
00283 if (!aval||!air||!ajc)
00284 { mexErrMsgTxt("DSDP: Problems with 2ND argument");}
00285 for (tnnz=0,jj=0;jj<nsubblocks;jj++){n=(int)blockn[jj];tnnz+=n*(n+1)/2;}
00286 if (nlhs>2){
00287 subs[0] = j; subs[1] = 0;
00288 index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
00289 X_cell_pr = mxCreateDoubleMatrix(tnnz,1,mxREAL);
00290 mxSetCell(X_OUT,index,X_cell_pr);
00291 xout=mxGetPr(X_cell_pr);
00292 if (tnnz>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
00293 }
00294
00295 for (ii=0; ii<=nvars; ii++){
00296 i=ii+1;
00297 if (i==nvars+1) i=0;
00298 tnnz=0; spot=ajc[ii]; blockj=sdpblockj;
00299 for (jj=0;jj<nsubblocks;jj++){
00300 n=(int)blockn[jj];
00301 if (sdpnmax<n) sdpnmax=n;
00302 if (ii==0){
00303 nn+=n;
00304 info=SDPConeSetBlockSize(sdpcone,blockj,n); CHKERR(info);
00305 info=SDPConeUsePackedFormat(sdpcone,blockj); CHKERR(info);
00306
00307 if (nlhs>2){info=SDPConeSetXArray(sdpcone,blockj,n,xout+tnnz,(n*n+n)/2); CHKERR(info);}
00308 info=CountNonzeroMatrices(blockn,nsubblocks,jj,nvars, air, ajc, &nzmats); CHKERR(info);
00309 info=SDPConeSetSparsity(sdpcone,blockj,nzmats); CHKERR(info);
00310 if (stat1<nzmats)stat1=nzmats;
00311 }
00312 for (tnnz2=tnnz+n*(n+1)/2,ijnnz=0;ijnnz<ajc[ii+1]-spot && air[spot+ijnnz]<tnnz2;ijnnz++){}
00313 if ( ijnnz==0 ){
00314 } else if(ijnnz==n*(n+1)/2){
00315 if (CheckForConstantMat(aval+spot,ijnnz,n)){
00316 info=SDPConeSetConstantMat(sdpcone,blockj,i,n,aval[spot]); CHKERR(info);
00317 } else {
00318 info=SDPConeSetADenseVecMat(sdpcone,blockj,i,n,1.0,aval+spot,ijnnz); CHKERR(info);
00319 }
00320 } else {
00321 info=SDPConeSetASparseVecMat(sdpcone,blockj,i,n,1.0,tnnz,air+spot,aval+spot,ijnnz); CHKERR(info);
00322 }
00323 tnnz+=n*(n+1)/2; spot+=ijnnz; blockj++;
00324 }
00325 }
00326 sdpblockj=sdpblockj+nsubblocks;
00327
00328 } else if (strcmp(conetype,"LB")==0 || strcmp(conetype,"UB")==0 ){
00329 subs[0] = j; subs[1] = 2;
00330 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00331 CA_cell_pr = mxGetCell(CA_IN,index);
00332 aval=mxGetPr(CA_cell_pr);
00333 info=DSDPCreateBCone(dsdp,&bcone); CHKERR(info);
00334 info=BConeAllocateBounds(bcone,nvars); CHKERR(info);
00335 for (i=0;i<nvars;i++){
00336 if (strcmp(conetype,"LB")==0){
00337 info=BConeSetLowerBound(bcone,i+1,aval[i]); CHKERR(info);
00338 } else {
00339 info=BConeSetUpperBound(bcone,i+1,aval[i]); CHKERR(info);
00340 }
00341 }
00342 if (nlhs>2){
00343 subs[0] = j; subs[1] = 0;
00344 index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
00345 X_cell_pr = mxCreateDoubleMatrix(nvars,1,mxREAL);
00346 mxSetCell(X_OUT,index,X_cell_pr);
00347 aval2=mxGetPr(X_cell_pr);
00348 info=BConeSetXArray(bcone,aval2,nvars); CHKERR(info);
00349 }
00350 } else if (strcmp(conetype,"LP")==0){
00351 subs[0] = j; subs[1] = 2;
00352 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00353 CA_cell_pr = mxGetCell(CA_IN,index);
00354 n = mxGetM(CA_cell_pr);
00355 if (lpnmax<n) lpnmax=n;
00356 nn+=n;
00357 aval=mxGetPr(CA_cell_pr); air =mxGetIr(CA_cell_pr); ajc =mxGetJc(CA_cell_pr);
00358 if (!aval||!air||!ajc)
00359 { mexErrMsgTxt("DSDP: Problems with 2ND argument");}
00360 info=DSDPCreateLPCone(dsdp,&lpcone); CHKERR(info);
00361 info=LPConeSetData2(lpcone,n,ajc,air,aval); CHKERR(info);
00362 if (nlhs>2){
00363 subs[0] = j; subs[1] = 0;
00364 index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
00365 X_cell_pr = mxCreateDoubleMatrix(n,1,mxREAL);
00366 mxSetCell(X_OUT,index,X_cell_pr);
00367 xout=mxGetPr(X_cell_pr);
00368 if (n>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
00369 info=LPConeSetXVec(lpcone,xout,n); CHKERR(info);
00370 }
00371
00372 } else if (strcmp(conetype,"FIXED")==0){
00373 int vari;
00374 subs[0] = j; subs[1] = 1;
00375 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00376 CA_cell_pr = mxGetCell(CA_IN,index);
00377 aval=mxGetPr(CA_cell_pr);
00378 it1 = mxGetM(CA_cell_pr);
00379 it2 = mxGetN(CA_cell_pr);
00380 subs[0] = j; subs[1] = 2;
00381 index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00382 CA_cell_pr = mxGetCell(CA_IN,index);
00383 aval2=mxGetPr(CA_cell_pr);
00384 if (nlhs>2){
00385 subs[0] = j; subs[1] = 0;
00386 index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
00387 X_cell_pr = mxCreateDoubleMatrix(it1*it2,1,mxREAL);
00388 mxSetCell(X_OUT,index,X_cell_pr);
00389 xout=mxGetPr(X_cell_pr);
00390 if (it1*it2>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
00391 } else {xout=0;}
00392 info=DSDPSetFixedVariables(dsdp,aval,aval2,xout,it1*it2); CHKERR(info);
00393 for (i=0;i<it1*it2;i++){
00394
00395
00396
00397
00398
00399 }
00400 }
00401 }
00402
00403
00404 if (y0){
00405 for (i=0;i<nvars;i++){ info = DSDPSetY0(dsdp,i+1,y0[i]); CHKERR(info);}
00406 }
00407
00408 reuse=(nvars-2)/sdpnmax;
00409 if (nvars<50 && reuse==0) reuse=1;
00410 if (reuse>=1) reuse++;
00411 reuse=reuse*reuse;
00412 if (nvars<2000 && reuse>10) reuse=10;
00413 if (reuse>12) reuse=12;
00414 info=DSDPReuseMatrix(dsdp,reuse); CHKERR(info);
00415
00416 info=DSDPGetPPObjective(dsdp,&zbar); CHKERR(info);
00417 info=DSDPGetR(dsdp,&r0); CHKERR(info);
00418 info=DSDPGetMaxIts(dsdp,&maxit); CHKERR(info);
00419 info=DSDPGetPenaltyParameter(dsdp,&penalty); CHKERR(info);
00420 info=DSDPGetPotentialParameter(dsdp,&rho); CHKERR(info);
00421 info=DSDPGetDualBound(dsdp,&dbound); CHKERR(info);
00422 info=DSDPGetGapTolerance(dsdp,&gaptol); CHKERR(info);
00423 info=DSDPGetRTolerance(dsdp,&inftol); CHKERR(info);
00424 info=DSDPGetBarrierParameter(dsdp,&mu0); CHKERR(info);
00425 info=DSDPGetMaxTrustRadius(dsdp,&maxtrust); CHKERR(info);
00426 info=DSDPGetStepTolerance(dsdp,&steptol); CHKERR(info);
00427 info=DSDPGetPTolerance(dsdp,&infptol); CHKERR(info);
00428 info=DSDPGetDataNorms(dsdp, datanorm); CHKERR(info);
00429 info=DSDPGetYBounds(dsdp,&ylow,&yhigh); CHKERR(info);
00430 if (datanorm[0]==0){info=DSDPSetYBounds(dsdp,-1.0,1.0);CHKERR(info);}
00431 if (nrhs >2){
00432 if(!mxIsStruct(PARS_IN)){
00433 mexErrMsgTxt("Fifth Parameter `OPTIONS' should be a structure.");}
00434
00435 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"maxit") ){
00436 maxit=(int) mxGetScalar(OPTIONS_FIELD);
00437 info=DSDPSetMaxIts(dsdp,maxit); CHKERR(info);}
00438 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"fastblas") ){
00439 fastblas= (int) mxGetScalar(OPTIONS_FIELD);}
00440 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"print") ){
00441 print_info= (int) mxGetScalar(OPTIONS_FIELD);
00442 printlevel=print_info;
00443 }
00444 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"bigM") ){
00445 rpos=(int)mxGetScalar(OPTIONS_FIELD);
00446 info=DSDPUsePenalty(dsdp,rpos); CHKERR(info);}
00447 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"penalty") ){
00448 penalty=mxGetScalar(OPTIONS_FIELD);
00449 info=DSDPSetPenaltyParameter(dsdp,penalty); CHKERR(info);}
00450 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"rho") ){
00451 rho= mxGetScalar(OPTIONS_FIELD);
00452 info=DSDPSetPotentialParameter(dsdp,rho); CHKERR(info);}
00453 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dynamicrho") ){
00454 drho= (int)mxGetScalar(OPTIONS_FIELD);
00455 info=DSDPUseDynamicRho(dsdp,drho); CHKERR(info);}
00456 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"zbar") ){
00457 zbar= mxGetScalar(OPTIONS_FIELD);
00458 info=DSDPSetZBar(dsdp,zbar); CHKERR(info);}
00459 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dual_bound") ){
00460 dbound= mxGetScalar(OPTIONS_FIELD);
00461 info=DSDPSetDualBound(dsdp,dbound); CHKERR(info);}
00462 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"reuse") ){
00463 reuse= (int)mxGetScalar(OPTIONS_FIELD);
00464 info=DSDPReuseMatrix(dsdp,reuse); CHKERR(info);}
00465 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"gaptol") ){
00466 gaptol= mxGetScalar(OPTIONS_FIELD);
00467 info=DSDPSetGapTolerance(dsdp,gaptol);CHKERR(info);}
00468 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"lp_barrier") ){
00469 lpb= mxGetScalar(OPTIONS_FIELD);
00470 if (lpb<0.1) lpb=0.1;
00471 if (lpcone){info=LPConeScaleBarrier(lpcone,lpb); CHKERR(info);}}
00472 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"lpb") ){
00473 lpb= mxGetScalar(OPTIONS_FIELD);
00474 if (lpb<0.1) lpb=0.1;
00475 if (lpcone){info=LPConeScaleBarrier(lpcone,lpb); CHKERR(info);}}
00476 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"cc") ){
00477 cc= mxGetScalar(OPTIONS_FIELD);
00478 info=DSDPAddObjectiveConstant(dsdp,cc); CHKERR(info);}
00479 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"inftol") ){
00480 inftol= mxGetScalar(OPTIONS_FIELD);
00481 info=DSDPSetRTolerance(dsdp,inftol); CHKERR(info);}
00482 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"infptol") ){
00483 infptol= mxGetScalar(OPTIONS_FIELD);
00484 info=DSDPSetPTolerance(dsdp,infptol);CHKERR(info);}
00485 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"pnormtol") ){
00486 pnormtol= mxGetScalar(OPTIONS_FIELD);
00487 info=DSDPSetPNormTolerance(dsdp,pnormtol);CHKERR(info);}
00488 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"boundy") ){
00489 yhigh= fabs(mxGetScalar(OPTIONS_FIELD)); ylow=-yhigh;
00490 info=DSDPSetYBounds(dsdp,ylow,yhigh);CHKERR(info);}
00491 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"r0") ){
00492 r0= mxGetScalar(OPTIONS_FIELD);
00493 info=DSDPSetR0(dsdp,r0); CHKERR(info);}
00494 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"mu0") ){
00495 mu0=mxGetScalar(OPTIONS_FIELD);
00496 info=DSDPSetBarrierParameter(dsdp,mu0);CHKERR(info);}
00497 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"maxtrustradius")){
00498 maxtrust= mxGetScalar(OPTIONS_FIELD);
00499 info=DSDPSetMaxTrustRadius(dsdp,maxtrust); CHKERR(info);}
00500 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"steptol") ){
00501 steptol = mxGetScalar(OPTIONS_FIELD);
00502 info=DSDPSetStepTolerance(dsdp,steptol); CHKERR(info);}
00503 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dobjmin") ){
00504 dtmp= mxGetScalar(OPTIONS_FIELD);
00505 info = DSDPSetDualLowerBound(dsdp,dtmp);}
00506 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dloginfo") ){
00507 iloginfo= (int) mxGetScalar(OPTIONS_FIELD);
00508 info=DSDPLogInfoAllow(iloginfo,0);}
00509 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"logtime") ){
00510 printsummary= (int) mxGetScalar(OPTIONS_FIELD);}
00511 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"printproblem") ){
00512 info=DSDPPrintData(dsdp,sdpcone,lpcone);}
00513 if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"xmaker") ){
00514 xmaker = (int) mxGetScalar(OPTIONS_FIELD);}
00515
00516
00517
00518
00519
00520
00521 }
00522
00523 info = DSDPSetMonitor(dsdp,DSDPPrintStats2,0); CHKERR(info);
00524
00525 info = DSDPSetup(dsdp);
00526 if (info){
00527 mexErrMsgTxt("DSDP: Setup Error, Probably out of memory");}
00528
00529
00530 info = DSDPSolve(dsdp); CHKERR(info);
00531 info = DSDPStopReason(dsdp,&reason);
00532 if (reason!=DSDP_INFEASIBLE_START){
00533 info=DSDPComputeX(dsdp);CHKERR(info);
00534 }
00535 info = DSDPGetSolutionType(dsdp,&pdfeasible); CHKERR(info);
00536
00537 if (printsummary){ DSDPEventLogSummary();}
00538
00539 if (info){
00540 mexErrMsgTxt("DSDP: Numerical error");}
00541
00542 if ( reason == DSDP_INFEASIBLE_START){
00543 mexErrMsgTxt("DSDP Terminated Due to Infeasible Starting Point\n");
00544 } else if (print_info){
00545
00546 if (reason == DSDP_CONVERGED)
00547 mexPrintf("DSDP Converged. \n");
00548 else if ( reason == DSDP_UPPERBOUND )
00549 mexPrintf("DSDP Converged: Dual Objective exceeds its bound\n");
00550 else if ( reason == DSDP_SMALL_STEPS )
00551 mexPrintf("DSDP Terminated Due to Small Steps\n");
00552 else if ( reason == DSDP_MAX_IT)
00553 mexPrintf("DSDP Terminated Due Maximum Number of Iterations\n");
00554 else if ( reason == DSDP_USER_TERMINATION)
00555 mexPrintf("DSDP Terminated By User\n");
00556 else if ( reason == DSDP_INFEASIBLE_START)
00557 mexPrintf("DSDP Terminated Due to Infeasible Starting Point\n");
00558 else
00559 mexPrintf("DSDP Finished.\n");
00560 }
00561
00562 if (pdfeasible==DSDP_UNBOUNDED){
00563 mexPrintf("DSDP: Dual Unbounded, Primal Infeasible\n");
00564 } else if (pdfeasible==DSDP_INFEASIBLE){
00565 mexPrintf("DSDP: Primal Unbounded, Dual Infeasible\n");
00566 }
00567
00568
00569 yout=mxGetPr(Y_OUT);
00570 info = DSDPGetY(dsdp,yout,nvars); CHKERR(info);
00571 if (info){
00572 mexErrMsgTxt("DSDP: Numerical error");}
00573
00574
00575
00576 if (STAT_OUT != NULL) mxDestroyArray(STAT_OUT) ;
00577 subs[0] = 1; subs[1] = 1;
00578 STAT_OUT = mxCreateStructArray(2,subs,nfields,fnames);
00579 info= DSDPGetDObjective(dsdp,&dobj); CHKERR(info);
00580 info= DSDPGetPObjective(dsdp,&pobj); CHKERR(info);
00581 info= DSDPGetR(dsdp,&dinf); CHKERR(info);
00582 info= DSDPStopReason(dsdp,&reason); CHKERR(info);
00583 info= DSDPGetIts(dsdp,&its); CHKERR(info);
00584 info= DSDPGetBarrierParameter(dsdp,&mu); CHKERR(info);
00585 info= DSDPGetStepLengths(dsdp,&pstep,&dstep); CHKERR(info);
00586 info= DSDPGetPnorm(dsdp,&pnorm); CHKERR(info);
00587 info= DSDPGetBarrierParameter(dsdp,&mu); CHKERR(info);
00588 info= DSDPGetYBounds(dsdp,&ylow,&yhigh); CHKERR(info);
00589
00590 if (pdfeasible==DSDP_UNBOUNDED){
00591 STAT_FIELD = mxCreateString("Unbounded");
00592 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
00593 } else if (pdfeasible==DSDP_INFEASIBLE){
00594 STAT_FIELD = mxCreateString("Infeasible");
00595 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
00596 } else {
00597 STAT_FIELD = mxCreateString("PDFeasible");
00598 mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
00599 }
00600
00601 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00602 stat=mxGetPr(STAT_FIELD); stat[0]=dobj;
00603 mxSetField(STAT_OUT,0,fnames[1],STAT_FIELD);
00604
00605 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00606 stat=mxGetPr(STAT_FIELD); stat[0]=pobj;
00607 mxSetField(STAT_OUT,0,fnames[2],STAT_FIELD);
00608
00609 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00610 stat=mxGetPr(STAT_FIELD); stat[0]=dobj;
00611 mxSetField(STAT_OUT,0,fnames[3],STAT_FIELD);
00612
00613 if (reason==DSDP_CONVERGED){
00614 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00615 stat=mxGetPr(STAT_FIELD); stat[0]=0;
00616 mxSetField(STAT_OUT,0,fnames[4],STAT_FIELD);
00617 } else {
00618 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00619 stat=mxGetPr(STAT_FIELD); stat[0]=(double)reason;
00620 if (stat[0]==0) stat[0]=-1;
00621 mxSetField(STAT_OUT,0,fnames[4],STAT_FIELD);
00622 }
00623
00624 if (reason==DSDP_CONVERGED){
00625 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00626 stat=mxGetPr(STAT_FIELD); stat[0]=0;
00627 if (pdfeasible==DSDP_UNBOUNDED){ stat[0]=1;}
00628 if (pdfeasible==DSDP_INFEASIBLE){ stat[0]=2;}
00629 mxSetField(STAT_OUT,0,fnames[5],STAT_FIELD);
00630 } else {
00631 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00632 stat=mxGetPr(STAT_FIELD); stat[0]=-1;
00633 if (reason==DSDP_INFEASIBLE_START)stat[0]=-6;
00634 if (reason==DSDP_UPPERBOUND)stat[0]=-5;
00635 if (reason==DSDP_MAX_IT)stat[0]=-3;
00636 mxSetField(STAT_OUT,0,fnames[5],STAT_FIELD);
00637 }
00638
00639 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00640 stat=mxGetPr(STAT_FIELD); stat[0]=(double) its;
00641 mxSetField(STAT_OUT,0,fnames[6],STAT_FIELD);
00642
00643 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00644 stat=mxGetPr(STAT_FIELD); stat[0]=dinf;
00645 mxSetField(STAT_OUT,0,fnames[7],STAT_FIELD);
00646
00647 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00648 stat=mxGetPr(STAT_FIELD); stat[0]=mu;
00649 mxSetField(STAT_OUT,0,fnames[8],STAT_FIELD);
00650
00651 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00652 stat=mxGetPr(STAT_FIELD); stat[0]=pstep;
00653 mxSetField(STAT_OUT,0,fnames[9],STAT_FIELD);
00654
00655 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00656 stat=mxGetPr(STAT_FIELD); stat[0]=dstep;
00657 mxSetField(STAT_OUT,0,fnames[10],STAT_FIELD);
00658
00659 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00660 stat=mxGetPr(STAT_FIELD); stat[0]=pnorm;
00661 mxSetField(STAT_OUT,0,fnames[11],STAT_FIELD);
00662
00663 itmp=100; if (its < itmp) itmp=its;
00664 STAT_FIELD = mxCreateDoubleMatrix(itmp, 1, mxREAL) ;
00665 stat=mxGetPr(STAT_FIELD);
00666 info= DSDPGetGapHistory(dsdp,stat,itmp); CHKERR(info);
00667 mxSetField(STAT_OUT,0,fnames[12],STAT_FIELD);
00668
00669 STAT_FIELD = mxCreateDoubleMatrix(itmp, 1, mxREAL) ;
00670 stat=mxGetPr(STAT_FIELD);
00671 info= DSDPGetRHistory(dsdp,stat,itmp); CHKERR(info);
00672 mxSetField(STAT_OUT,0,fnames[13],STAT_FIELD);
00673
00674 STAT_FIELD = mxCreateDoubleMatrix(6, 1, mxREAL) ;
00675 stat=mxGetPr(STAT_FIELD);
00676 info = DSDPGetFinalErrors(dsdp,stat); CHKERR(info);
00677 mxSetField(STAT_OUT,0,fnames[14],STAT_FIELD);
00678
00679 STAT_FIELD = mxCreateDoubleMatrix(3, 1, mxREAL) ;
00680 stat=mxGetPr(STAT_FIELD);
00681 info = DSDPGetDataNorms(dsdp,stat); CHKERR(info);
00682 dtmp=stat[0];stat[0]=stat[1];stat[1]=dtmp;
00683 mxSetField(STAT_OUT,0,fnames[15],STAT_FIELD);
00684
00685 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00686 stat=mxGetPr(STAT_FIELD);
00687 info = DSDPGetYMaxNorm(dsdp,stat); CHKERR(info);
00688 mxSetField(STAT_OUT,0,fnames[16],STAT_FIELD);
00689
00690 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00691 stat=mxGetPr(STAT_FIELD);
00692 stat[0]=yhigh;
00693 mxSetField(STAT_OUT,0,fnames[17],STAT_FIELD);
00694
00695 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00696 stat=mxGetPr(STAT_FIELD);
00697 info = DSDPGetPenaltyParameter(dsdp,stat); CHKERR(info);
00698 mxSetField(STAT_OUT,0,fnames[18],STAT_FIELD);
00699
00700 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00701 stat=mxGetPr(STAT_FIELD);
00702 info = DSDPGetTraceX(dsdp,stat); CHKERR(info);
00703 mxSetField(STAT_OUT,0,fnames[19],STAT_FIELD);
00704
00705 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00706 info = DSDPGetReuseMatrix(dsdp,&its); CHKERR(info);
00707 stat=mxGetPr(STAT_FIELD); stat[0]=(double) its;
00708 mxSetField(STAT_OUT,0,fnames[20],STAT_FIELD);
00709
00710 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00711 stat=mxGetPr(STAT_FIELD);
00712 info = DSDPGetPotentialParameter(dsdp,stat); CHKERR(info);
00713 mxSetField(STAT_OUT,0,fnames[21],STAT_FIELD);
00714
00715 if (xmaker){
00716 STAT_FIELD = mxCreateDoubleMatrix(nvars+1, 1, mxREAL) ;
00717 stat=mxGetPr(STAT_FIELD);
00718 info = DSDPGetYMakeX(dsdp,stat,nvars+1); CHKERR(info);
00719 mxSetField(STAT_OUT,0,fnames[22],STAT_FIELD);
00720
00721 STAT_FIELD = mxCreateDoubleMatrix(nvars+1, 1, mxREAL) ;
00722 stat=mxGetPr(STAT_FIELD);
00723 info = DSDPGetDYMakeX(dsdp,stat,nvars+1); CHKERR(info);
00724 mxSetField(STAT_OUT,0,fnames[23],STAT_FIELD);
00725
00726 STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00727 stat=mxGetPr(STAT_FIELD);
00728 info = DSDPGetMuMakeX(dsdp,stat); CHKERR(info);
00729 mxSetField(STAT_OUT,0,fnames[24],STAT_FIELD);
00730 }
00731
00732
00733 info = DSDPDestroy(dsdp); CHKERR(info);
00734
00735 return;
00736 }
00737
00738
00739
00740 #undef __FUNCT__
00741 #define __FUNCT__ "CheckForConstantMat"
00742 static int CheckForConstantMat(double*v,int nnz, int n){
00743 int i;double vv;
00744 if (n<=1){ return 0; }
00745 if (nnz!=(n*n+n)/2){ return 0; }
00746 for (vv=v[0],i=1;i<nnz;i++){
00747 if (v[i]!=vv){ return 0;}
00748 }
00749 return 1;
00750 }
00751
00752 #undef __FUNCT__
00753 #define __FUNCT__ "CountNonzeroMatrices"
00754 static int CountNonzeroMatrices(double *blocksize, int nblocks, int block, int nvars, int *indd, int*nnnz, int *nnzmats){
00755 int i,j,n,nzmats=0;
00756 int marker1=0,marker2=0;
00757 for (i=0;i<block;i++){n=(int)blocksize[i]; marker1=marker1+n*(n+1)/2;}
00758 n=(int)blocksize[block];
00759 marker2=marker1+n*(n+1)/2;
00760 for (i=0;i<nvars;i++){
00761 j=nnnz[i];
00762 while (indd[j]<marker1 && j< nnnz[i+1]){j++;}
00763 if (j<nnnz[i+1] && indd[j]<marker2){ nzmats++;}
00764 }
00765 *nnzmats=nzmats;
00766 return 0;
00767 }
00768
00769
00770 static int DSDPPrintStats2(DSDP dsdp, void* dummy){
00771
00772 int iter,info;
00773 double pobj,dobj,pstp=0,dstp,mu,res,pnorm,pinfeas;
00774 DSDPTerminationReason reason;
00775
00776 if(printlevel<=0) return(0);
00777
00778 info = DSDPStopReason(dsdp,&reason);
00779 info = DSDPGetIts(dsdp,&iter);
00780
00781 if( (reason!=CONTINUE_ITERATING) || ((iter % printlevel)==0)){
00782
00783 info = DSDPGetDDObjective(dsdp,&dobj);
00784 info = DSDPGetPPObjective(dsdp,&pobj);
00785 info = DSDPGetR(dsdp,&res);
00786 info = DSDPGetPInfeasibility(dsdp,&pinfeas);
00787 info = DSDPGetStepLengths(dsdp,&pstp,&dstp);
00788 info = DSDPGetBarrierParameter(dsdp,&mu);
00789 info = DSDPGetPnorm(dsdp,&pnorm);
00790
00791
00792 if (iter==0){
00793 mexPrintf("Iter PP Objective DD Objective PInfeas DInfeas Nu StepLength Pnrm\n")
00794 ;
00795 mexPrintf("----------------------------------------------------------------------------------------\n")
00796 ;
00797 }
00798 mexPrintf("%-3d %16.8e %16.8e %9.1e %9.1e %9.1e",iter,pobj,dobj,pinfeas,res,mu);
00799 mexPrintf(" %4.2f %4.2f",pstp,dstp);
00800 if (pnorm>1.0e3){
00801 mexPrintf(" %1.0e \n",pnorm);
00802 } else {
00803 mexPrintf(" %5.2f \n",pnorm);
00804 }
00805 fflush(NULL);
00806 }
00807 return(0);
00808 }