00001 #include "dsdp5.h"
00002 #include <stdio.h>
00003 #include <string.h>
00004 #include <math.h>
00005 #include <stdlib.h>
00006
00011 char help[]="\n\
00012 Compute the Lovasz theta number for a graph. This number is an upper bound for \n\
00013 the maximum clique of a graph, a lower bound for the mimimal graph coloring, and serves \n\
00014 as a bound for several other combitorial graph problems. The number is the solution to \n\
00015 a semidfinite program. \n\n\
00016 The file should be the complement of the graph \n\
00017 This file also demonstrates the use of customized data matrices in DSDP.\n\n\
00018 DSDP Usage: theta <graph complement filename> \n";
00019
00020 typedef struct {
00021 int v1,v2;
00022 } EdgeMat;
00023
00024 #include "../src/sdp/dsdpdatamat_impl.h"
00025 extern int SDPConeAddDataMatrix(SDPCone,int, int, int, char, struct DSDPDataMat_Ops*, void*);
00026
00027
00028
00029 static struct DSDPDataMat_Ops edgematop;
00030 static int EdgeMatOperationsInitialize(struct DSDPDataMat_Ops*);
00031
00032 static struct DSDPDataMat_Ops onematops;
00033 static int OneMatOpsInitialize(struct DSDPDataMat_Ops*);
00034
00035 static struct DSDPDataMat_Ops identitymatops;
00036 static int IdentitymatOperationsInitialize(struct DSDPDataMat_Ops*);
00037
00038 static int ReadGraphFromFile(char*,int*, int*, EdgeMat*[]);
00039 int SetThetaData(DSDP, SDPCone, int, int, EdgeMat[]);
00040 int LovaszTheta(int argc,char *argv[]);
00041
00042 int main(int argc,char *argv[]){
00043 int info;
00044 info=LovaszTheta(argc,argv);
00045 return 0;
00046 }
00047
00056 int LovaszTheta(int argc,char *argv[]){
00057
00058 int info,kk,nedges,nnodes;
00059 EdgeMat *Edges;
00060 DSDP dsdp;
00061 SDPCone sdpcone;
00062
00063 if (argc<2){ printf("%s",help); return(1); }
00064
00065 info = ReadGraphFromFile(argv[1],&nnodes,&nedges,&Edges);
00066 if (info){ printf("Problem reading file\n"); return 1; }
00067
00068 info = DSDPCreate(nedges+1,&dsdp);
00069 info = DSDPCreateSDPCone(dsdp,1,&sdpcone);
00070 info = SDPConeSetBlockSize(sdpcone,0,nnodes);
00071 info = SDPConeUsePackedFormat(sdpcone,0);
00072 info = SDPConeSetSparsity(sdpcone,0,nedges+1);
00073 if (info){ printf("Out of memory\n"); return 1; }
00074 info = SetThetaData(dsdp, sdpcone, nnodes, nedges, Edges);
00075 if (info){ printf("Out of memory\n"); return 1; }
00076
00077 info = DSDPSetGapTolerance(dsdp,0.001);
00078 info = DSDPSetZBar(dsdp,nnodes+1);
00079 info = DSDPReuseMatrix(dsdp,10);
00080
00081 for (kk=1; kk<argc-1; kk++){
00082 if (strncmp(argv[kk],"-dloginfo",8)==0){
00083 info=DSDPLogInfoAllow(DSDP_TRUE,0);
00084 } else if (strncmp(argv[kk],"-params",7)==0){
00085 info=DSDPReadOptions(dsdp,argv[kk+1]);
00086 } else if (strncmp(argv[kk],"-help",5)==0){
00087 printf("%s\n",help);
00088 }
00089 }
00090 info=DSDPSetOptions(dsdp,argv,argc);
00091
00092 if (info){ printf("Out of memory\n"); return 1; }
00093 info = DSDPSetStandardMonitor(dsdp,1);
00094 if (info){ printf("Monitor Problem \n"); return 1; }
00095
00096 info = DSDPSetup(dsdp);
00097 if (info){ printf("Out of memory\n"); return 1; }
00098 if (0==1){info=SDPConeCheckData(sdpcone);}
00099
00100 info = DSDPSolve(dsdp);
00101 if (info){ printf("Numberical error\n"); return 1; }
00102 info=DSDPComputeX(dsdp);DSDPCHKERR(info);
00103
00104 if (0==1){
00105 int n; double *xx;
00106 info=SDPConeGetXArray(sdpcone,0,&xx,&n);
00107 }
00108
00109 info = DSDPDestroy(dsdp);
00110 free(Edges);
00111
00112 return 0;
00113 }
00114
00126 int SetThetaData(DSDP dsdp, SDPCone sdpcone, int nodes, int edges, EdgeMat Edge[]){
00127
00128 int i,info;
00129
00130
00131
00132
00133 info=OneMatOpsInitialize(&onematops);
00134 info=SDPConeAddDataMatrix(sdpcone,0,0,nodes,'P',&onematops,0);
00135
00136
00137 info=EdgeMatOperationsInitialize(&edgematop);
00138 for (i=0; i<edges; i++){
00139 info = SDPConeAddDataMatrix(sdpcone,0,i+1,nodes,'P',&edgematop,(void*)&Edge[i]);
00140 info = DSDPSetDualObjective(dsdp,i+1,0.0);
00141 info = DSDPSetY0(dsdp,i+1,0.0);
00142 if (info) return info;
00143 }
00144
00145
00146 info = IdentitymatOperationsInitialize(&identitymatops);
00147 info = SDPConeAddDataMatrix(sdpcone,0,edges+1,nodes,'P',&identitymatops,0);
00148 info = DSDPSetDualObjective(dsdp,edges+1,1.0);
00149 info = DSDPSetY0(dsdp,edges+1,-10*nodes-1);
00150
00151
00152 info = DSDPSetR0(dsdp,0);
00153
00154 return(0);
00155 }
00156
00157 #define BUFFERSIZ 100
00158
00159 #undef __FUNCT__
00160 #define __FUNCT__ "ParseGraphline"
00161 static int ParseGraphline(char * thisline,int *row,int *col,double *value,
00162 int *gotem){
00163
00164 int temp;
00165 int rtmp,coltmp;
00166
00167 rtmp=-1, coltmp=-1, *value=0.0;
00168 temp=sscanf(thisline,"%d %d %lf",&rtmp,&coltmp,value);
00169 if (temp==3 && rtmp>0 && coltmp>0) *gotem=3;
00170 else if (temp==2 && rtmp>0 && coltmp>0){ *value = 1.0; *gotem=3;}
00171 else *gotem=0;
00172 *row=rtmp-1; *col=coltmp-1;
00173 if (*gotem && (*col < 0 || *row < 0)){
00174 printf("Node Number must be positive.\n, %s\n",thisline);
00175 }
00176 return(0);
00177 }
00178
00179 #undef __FUNCT__
00180 #define __FUNCT__ "ReadGraphFromFile"
00181 static int ReadGraphFromFile(char* filename,int *nnodes, int *nedges, EdgeMat**EE){
00182
00183 FILE*fp;
00184 char thisline[BUFFERSIZ]="*";
00185 int i,k=0,line=0,nodes,edges,gotone=3;
00186 int info,row,col;
00187 double value;
00188 EdgeMat *E;
00189
00190 fp=fopen(filename,"r");
00191 if (!fp){ printf("Cannot open file %s !",filename); return(1); }
00192
00193 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){
00194 fgets(thisline,BUFFERSIZ,fp); line++; }
00195
00196 if (sscanf(thisline,"%d %d",&nodes, &edges)!=2){
00197 printf("First line must contain the number of nodes and number of edges\n");
00198 return 1;
00199 }
00200
00201 E=(EdgeMat*)malloc(edges*sizeof(EdgeMat)); *EE=E;
00202 for (i=0; i<edges; i++){ E[i].v1=0; E[i].v2=0; }
00203
00204 while(!feof(fp) && gotone){
00205 thisline[0]='\0';
00206 fgets(thisline,BUFFERSIZ,fp); line++;
00207 info = ParseGraphline(thisline,&row,&col,&value,&gotone);
00208 if (gotone && k<edges &&
00209 col < nodes && row < nodes && col >= 0 && row >= 0){
00210 if (row > col){i=row;row=col;col=i;}
00211 if (row == col){}
00212 else { E[k].v1=row; E[k].v2=col; k++;}
00213 } else if (gotone && k>=edges) {
00214 printf("To many edges in file.\nLine %d, %s\n",line,thisline);
00215 return 1;
00216 } else if (gotone&&(col >= nodes || row >= nodes || col < 0 || row < 0)){
00217 printf("Invalid node number.\nLine %d, %s\n",line,thisline);
00218 return 1;
00219 }
00220 }
00221
00222 *nnodes=nodes; *nedges=k;
00223
00224 return 0;
00225 }
00226
00227
00228
00229 static int EdgeMatDestroy(void*);
00230 static int EdgeMatView(void*);
00231 static int EdgeMatVecVec(void*, double[], int, double *);
00232 static int EdgeMatDot(void*, double[], int, int, double *);
00233 static int EdgeMatGetRank(void*, int*, int);
00234 static int EdgeMatFactor(void*);
00235 static int EdgeMatGetEig(void*, int, double*, double[], int, int[], int*);
00236 static int EdgeMatAddRowMultiple(void*, int, double, double[], int);
00237 static int EdgeMatAddMultiple(void*, double, double[], int, int);
00238 static int EdgeMatGetRowNnz(void*, int, int[], int*, int);
00239
00240 static int EdgeMatDestroy(void* AA){
00241 return 0;
00242 }
00243 static int EdgeMatVecVec(void* A, double x[], int n, double *v){
00244 EdgeMat*E=(EdgeMat*)A;
00245 *v=2*x[E->v1]*x[E->v2];
00246 return 0;
00247 }
00248 static int EdgeMatDot(void* A, double x[], int nn, int n, double *v){
00249 EdgeMat*E=(EdgeMat*)A;
00250 int k=E->v2*(E->v2+1)/2 + E->v1;
00251 *v=2*x[k];
00252 return 0;
00253 }
00254 static int EdgeMatView(void* A){
00255 EdgeMat*E=(EdgeMat*)A;
00256 printf(" Row: %d, Column: %d\n",E->v1,E->v2);
00257 printf(" Row: %d, Column: %d\n",E->v2,E->v1);
00258 return 0;
00259 }
00260 static int EdgeMatFactor(void* A){
00261 return 0;
00262 }
00263 static int EdgeMatGetRank(void *A, int*rank, int n){
00264 *rank=2;
00265 return 0;
00266 }
00267 static int EdgeMatGetEig(void*A, int neig, double *eig, double v[], int n,int spind[], int *nind){
00268 EdgeMat*E=(EdgeMat*)A;
00269 double tt=1.0/sqrt(2.0);
00270 memset((void*)v,0,(n)*sizeof(double));
00271 memset((void*)spind,0,(n)*sizeof(int));
00272 if (neig==0){
00273 v[E->v1]=tt;v[E->v2]=tt;*eig=1;
00274 spind[0]=E->v1;spind[1]=E->v2; *nind=2;
00275 } else if (neig==1){
00276 v[E->v1]=-tt;v[E->v2]=tt;*eig=-1;
00277 spind[0]=E->v1;spind[1]=E->v2; *nind=2;
00278 } else { *eig=0;*nind=0;}
00279 return 0;
00280 }
00281 static int EdgeMatGetRowNnz(void*A, int nrow, int nz[], int *nnzz, int n){
00282 EdgeMat*E=(EdgeMat*)A;
00283 if (nrow==E->v1){ nz[E->v2]++; *nnzz=1;}
00284 else if (nrow==E->v2){nz[E->v1]++; *nnzz=1;}
00285 else {*nnzz=0;}
00286 return 0;
00287 }
00288 static int EdgeMatAddRowMultiple(void*A, int nrow, double dd, double rrow[], int n){
00289 EdgeMat*E=(EdgeMat*)A;
00290 if (nrow==E->v1){ rrow[E->v2]+=dd;}
00291 else if (nrow==E->v2){rrow[E->v1]+=dd;}
00292 return 0;
00293 }
00294 static int EdgeMatAddMultiple(void*A, double dd, double vv[], int nn, int n){
00295 EdgeMat*E=(EdgeMat*)A;
00296 int k=E->v2*(E->v2+1)/2 + E->v1;
00297 vv[k]+=dd;
00298 return 0;
00299 }
00300 static int EdgeMatFNorm(void*A, int n, double *fnorm){
00301 *fnorm=2.0;
00302 return 0;
00303 }
00304 static int EdgeMatCountNonzeros(void*A, int *nnz, int n){
00305 *nnz=1;
00306 return 0;
00307 }
00308 static const char *datamatname="THETA EDGE MATRIX";
00309 static int EdgeMatOperationsInitialize(struct DSDPDataMat_Ops* edgematoperator){
00310 int info;
00311 if (edgematoperator==NULL) return 0;
00312 info=DSDPDataMatOpsInitialize(edgematoperator); if (info){ return info;}
00313 edgematoperator->matfactor1=EdgeMatFactor;
00314 edgematoperator->matgetrank=EdgeMatGetRank;
00315 edgematoperator->matgeteig=EdgeMatGetEig;
00316 edgematoperator->matvecvec=EdgeMatVecVec;
00317 edgematoperator->matrownz=EdgeMatGetRowNnz;
00318 edgematoperator->matdot=EdgeMatDot;
00319 edgematoperator->matfnorm2=EdgeMatFNorm;
00320 edgematoperator->matnnz=EdgeMatCountNonzeros;
00321 edgematoperator->mataddrowmultiple=EdgeMatAddRowMultiple;
00322 edgematoperator->mataddallmultiple=EdgeMatAddMultiple;
00323 edgematoperator->matdestroy=EdgeMatDestroy;
00324 edgematoperator->matview=EdgeMatView;
00325 edgematoperator->matname=datamatname;
00326 edgematoperator->id=25;
00327 return 0;
00328 }
00329
00330
00331 static int OneMatDestroy(void*);
00332 static int OneMatView(void*);
00333 static int OneMatVecVec(void*, double[], int, double *);
00334 static int OneMatDot(void*, double[], int, int, double *);
00335 static int OneMatGetRank(void*, int*, int);
00336 static int OneMatFactor(void*);
00337 static int OneMatGetEig(void*, int, double*, double[], int, int[], int*);
00338 static int OneMatRowNnz(void*, int, int[], int*, int);
00339 static int OneMatAddRowMultiple(void*, int, double, double[], int);
00340 static int OneMatAddMultiple(void*, double, double[], int,int);
00341
00342
00343 static int OneMatFactor(void*A){return 0;}
00344 static int OneMatGetRank(void *A, int *rank, int n){*rank=1;return 0;}
00345 static int OneMatFNorm2(void*AA, int n, double *fnorm2){*fnorm2=1.0*n*n;return 0;}
00346 static int OneMatCountNonzeros(void*AA, int *nnz, int n){*nnz=n*n;return 0;}
00347 static int OneMatDot(void* A, double x[], int nn, int n, double *v){
00348 double dtmp=0.0;
00349 int i,j;
00350 for (i=0;i<n;i++){
00351 for (j=0;j<i;j++,x++){dtmp+= (*x);}
00352 dtmp+= (*x);x++;
00353 }
00354 *v=-2*dtmp;
00355 return 0;
00356 }
00357 static int OneMatVecVec(void* A, double x[], int n, double *v){
00358 double dtmp=0.0;
00359 int i;
00360 for (i=0; i<n; i++){dtmp+=x[i];}
00361 *v=-dtmp*dtmp;
00362 return 0;
00363 }
00364 static int OneMatAddMultiple(void*A, double ddd, double vv[], int nn, int n){
00365 int i,j;
00366 for (i=0;i<n;i++){
00367 for (j=0;j<i;j++,vv++){(*vv)+=-ddd;}
00368 (*vv)+=-ddd; vv++;
00369 }
00370 return 0;
00371 }
00372 static int OneMatAddRowMultiple(void*A, int nrow, double ddd, double row[], int n){
00373 int i;
00374 for (i=0;i<n;i++){row[i] -= ddd;}
00375 row[nrow] -= ddd;
00376 return 0;
00377 }
00378 static int OneMatGetEig(void*A, int neig, double *eig, double v[], int n, int spind[], int *nind){
00379 int i;
00380 if (neig==0){ *eig=-1; for (i=0;i<n;i++){ v[i]=1.0; spind[i]=i;} *nind=n;
00381 } else { *eig=0; for (i=0;i<n;i++){ v[i]=0.0; } *nind=0;
00382 }
00383 return 0;
00384 }
00385 static int OneMatRowNnz(void*A, int row, int nz[], int *nnz, int n){
00386 int i;
00387 for (i=0;i<n;i++){ nz[i]++; }
00388 *nnz=n;
00389 return 0;
00390 }
00391 static int OneMatView(void* AA){
00392 printf("Every element of the matrix is the same: -1\n");
00393 return 0;
00394 }
00395 static int OneMatDestroy(void* A){return 0;}
00396
00397 static const char *mat1name="THETA ALL ELEMENTS EQUAL -ONE";
00398 static int OneMatOpsInitialize(struct DSDPDataMat_Ops* mat1ops){
00399 int info;
00400 if (mat1ops==NULL) return 0;
00401 info=DSDPDataMatOpsInitialize(mat1ops); DSDPCHKERR(info);
00402 mat1ops->matfactor1=OneMatFactor;
00403 mat1ops->matgetrank=OneMatGetRank;
00404 mat1ops->matgeteig=OneMatGetEig;
00405 mat1ops->matvecvec=OneMatVecVec;
00406 mat1ops->matdot=OneMatDot;
00407 mat1ops->mataddrowmultiple=OneMatAddRowMultiple;
00408 mat1ops->mataddallmultiple=OneMatAddMultiple;
00409 mat1ops->matdestroy=OneMatDestroy;
00410 mat1ops->matview=OneMatView;
00411 mat1ops->matrownz=OneMatRowNnz;
00412 mat1ops->matfnorm2=OneMatFNorm2;
00413 mat1ops->matnnz=OneMatCountNonzeros;
00414 mat1ops->id=18;
00415 mat1ops->matname=mat1name;
00416 return 0;
00417 }
00418
00419
00420 static int IdentityMatDestroy(void*);
00421 static int IdentityMatView(void*);
00422 static int IdentityMatVecVec(void*, double[], int, double *);
00423 static int IdentityMatDot(void*, double[], int, int, double *);
00424 static int IdentityMatGetRank(void*, int*,int);
00425 static int IdentityMatFactor(void*);
00426 static int IdentityMatGetEig(void*, int, double*, double[], int, int[], int*);
00427 static int IdentityMatAddRowMultiple(void*, int, double, double[], int);
00428 static int IdentityMatAddMultiple(void*, double, double[], int, int);
00429 static int IdentityMatGetRowNnz(void*, int, int[], int*, int);
00430
00431 static int IdentityMatDestroy(void* AA){return 0;}
00432 static int IdentityMatFNorm2(void* AA, int n, double *v){*v=1.0*n;return 0;}
00433 static int IdentityMatGetRank(void *AA, int*rank, int n){*rank=n;return 0;}
00434 static int IdentityMatFactor(void*A){return 0;}
00435 static int IdentityMatCountNonzeros(void*A, int *nnz, int n){*nnz=n;return 0;}
00436 static int IdentityMatVecVec(void* AA, double x[], int n, double *v){
00437 int i;
00438 *v=0;
00439 for (i=0;i<n;i++){ *v+=x[i]*x[i]; }
00440 return 0;
00441 }
00442 static int IdentityMatDot(void* AA, double x[], int nn, int n, double *v){
00443 int i;
00444 double vv=0;
00445 for (i=0;i<n;i++){ vv+=x[((i+1)*(i+2))/2-1];}
00446 *v = 2*vv;
00447 return 0;
00448 }
00449 static int IdentityMatView(void* AA){
00450 printf("Identity matrix: All Diagonal elements equal 1.0\n");
00451 return 0;
00452 }
00453 static int IdentityMatGetEig(void*AA, int neig, double *eig, double v[], int n, int spind[], int *nind){
00454 if (neig<0 || neig>=n){ *eig=0; return 0;}
00455 memset((void*)v,0,(n)*sizeof(double));
00456 *eig=1.0; v[neig]=1.0; spind[0]=neig; *nind=1;
00457 return 0;
00458 }
00459 static int IdentityMatGetRowNnz(void*A, int nrow, int nz[], int *nnzz, int n){
00460 if (nrow>=0 && nrow < n){
00461 *nnzz=1; nz[nrow]++;
00462 } else { *nnzz=0;
00463 }
00464 return 0;
00465 }
00466 static int IdentityMatAddRowMultiple(void*A, int nrow, double dd, double rrow[], int n){
00467 rrow[nrow] += dd;return 0;
00468 }
00469 static int IdentityMatAddMultiple(void*A, double dd, double vv[], int nn, int n){
00470 int i;
00471 for (i=0;i<n;i++){ vv[(i+1)*(i+2)/2-1] += dd;}
00472 return 0;
00473 }
00474
00475 static const char *eyematname="THETA IDENTITY MATRIX";
00476 static int IdentitymatOperationsInitialize(struct DSDPDataMat_Ops* imatops){
00477 int info;
00478 if (imatops==NULL) return 0;
00479 info=DSDPDataMatOpsInitialize(imatops); if (info){return info;}
00480 imatops->matfactor1=IdentityMatFactor;
00481 imatops->matgetrank=IdentityMatGetRank;
00482 imatops->matgeteig=IdentityMatGetEig;
00483 imatops->matvecvec=IdentityMatVecVec;
00484 imatops->matrownz=IdentityMatGetRowNnz;
00485 imatops->matdot=IdentityMatDot;
00486 imatops->matfnorm2=IdentityMatFNorm2;
00487 imatops->matnnz=IdentityMatCountNonzeros;
00488 imatops->mataddrowmultiple=IdentityMatAddRowMultiple;
00489 imatops->mataddallmultiple=IdentityMatAddMultiple;
00490 imatops->matdestroy=IdentityMatDestroy;
00491 imatops->matview=IdentityMatView;
00492 imatops->id=12;
00493 imatops->matname=eyematname;
00494 return 0;
00495 }