12Compute the Lovasz theta number for a graph. This number is an upper bound for \n\
13the maximum clique of a graph, a lower bound for the mimimal graph coloring, and serves \n\
14as a bound for several other combitorial graph problems. The number is the solution to \n\
15a semidfinite program. \n\n\
16The file should be the complement of the graph \n\
17This file also demonstrates the use of customized data matrices in DSDP.\n\n\
18DSDP Usage: theta <graph complement filename> \n";
24#include "../src/sdp/dsdpdatamat_impl.h"
38static int ReadGraphFromFile(
char*,
int*,
int*, EdgeMat*[]);
39int SetThetaData(
DSDP,
SDPCone,
int,
int, EdgeMat[]);
40int LovaszTheta(
int argc,
char *argv[]);
42int main(
int argc,
char *argv[]){
44 info=LovaszTheta(argc,argv);
56int LovaszTheta(
int argc,
char *argv[]){
58 int info,kk,nedges,nnodes;
63 if (argc<2){ printf(
"%s",help);
return(1); }
65 info = ReadGraphFromFile(argv[1],&nnodes,&nedges,&Edges);
66 if (info){ printf(
"Problem reading file\n");
return 1; }
68 info = DSDPCreate(nedges+1,&dsdp);
69 info = DSDPCreateSDPCone(dsdp,1,&sdpcone);
70 info = SDPConeSetBlockSize(sdpcone,0,nnodes);
71 info = SDPConeUsePackedFormat(sdpcone,0);
72 info = SDPConeSetSparsity(sdpcone,0,nedges+1);
73 if (info){ printf(
"Out of memory\n");
return 1; }
74 info = SetThetaData(dsdp, sdpcone, nnodes, nedges, Edges);
75 if (info){ printf(
"Out of memory\n");
return 1; }
77 info = DSDPSetGapTolerance(dsdp,0.001);
78 info = DSDPSetZBar(dsdp,nnodes+1);
79 info = DSDPReuseMatrix(dsdp,10);
81 for (kk=1; kk<argc-1; kk++){
82 if (strncmp(argv[kk],
"-dloginfo",8)==0){
84 }
else if (strncmp(argv[kk],
"-params",7)==0){
85 info=DSDPReadOptions(dsdp,argv[kk+1]);
86 }
else if (strncmp(argv[kk],
"-help",5)==0){
90 info=DSDPSetOptions(dsdp,argv,argc);
92 if (info){ printf(
"Out of memory\n");
return 1; }
93 info = DSDPSetStandardMonitor(dsdp,1);
94 if (info){ printf(
"Monitor Problem \n");
return 1; }
96 info = DSDPSetup(dsdp);
97 if (info){ printf(
"Out of memory\n");
return 1; }
98 if (0==1){info=SDPConeCheckData(sdpcone);}
100 info = DSDPSolve(dsdp);
101 if (info){ printf(
"Numberical error\n");
return 1; }
102 info=DSDPComputeX(dsdp);DSDPCHKERR(info);
106 info=SDPConeGetXArray(sdpcone,0,&xx,&n);
109 info = DSDPDestroy(dsdp);
126int SetThetaData(
DSDP dsdp,
SDPCone sdpcone,
int nodes,
int edges, EdgeMat Edge[]){
133 info=OneMatOpsInitialize(&onematops);
134 info=SDPConeAddDataMatrix(sdpcone,0,0,nodes,
'P',&onematops,0);
137 info=EdgeMatOperationsInitialize(&edgematop);
138 for (i=0; i<edges; i++){
139 info = SDPConeAddDataMatrix(sdpcone,0,i+1,nodes,
'P',&edgematop,(
void*)&Edge[i]);
140 info = DSDPSetDualObjective(dsdp,i+1,0.0);
141 info = DSDPSetY0(dsdp,i+1,0.0);
142 if (info)
return info;
146 info = IdentitymatOperationsInitialize(&identitymatops);
147 info = SDPConeAddDataMatrix(sdpcone,0,edges+1,nodes,
'P',&identitymatops,0);
148 info = DSDPSetDualObjective(dsdp,edges+1,1.0);
149 info = DSDPSetY0(dsdp,edges+1,-10*nodes-1);
152 info = DSDPSetR0(dsdp,0);
160#define __FUNCT__ "ParseGraphline"
161static int ParseGraphline(
char * thisline,
int *row,
int *col,
double *value,
167 rtmp=-1, coltmp=-1, *value=0.0;
168 temp=sscanf(thisline,
"%d %d %lf",&rtmp,&coltmp,value);
169 if (temp==3 && rtmp>0 && coltmp>0) *gotem=3;
170 else if (temp==2 && rtmp>0 && coltmp>0){ *value = 1.0; *gotem=3;}
172 *row=rtmp-1; *col=coltmp-1;
173 if (*gotem && (*col < 0 || *row < 0)){
174 printf(
"Node Number must be positive.\n, %s\n",thisline);
180#define __FUNCT__ "ReadGraphFromFile"
181static int ReadGraphFromFile(
char* filename,
int *nnodes,
int *nedges, EdgeMat**EE){
184 char thisline[BUFFERSIZ]=
"*";
185 int i,k=0,line=0,nodes,edges,gotone=3;
190 fp=fopen(filename,
"r");
191 if (!fp){ printf(
"Cannot open file %s !",filename);
return(1); }
193 while(!feof(fp) && (thisline[0] ==
'*' || thisline[0] ==
'"') ){
194 fgets(thisline,BUFFERSIZ,fp); line++; }
196 if (sscanf(thisline,
"%d %d",&nodes, &edges)!=2){
197 printf(
"First line must contain the number of nodes and number of edges\n");
201 E=(EdgeMat*)malloc(edges*
sizeof(EdgeMat)); *EE=E;
202 for (i=0; i<edges; i++){ E[i].v1=0; E[i].v2=0; }
204 while(!feof(fp) && gotone){
206 fgets(thisline,BUFFERSIZ,fp); line++;
207 info = ParseGraphline(thisline,&row,&col,&value,&gotone);
208 if (gotone && k<edges &&
209 col < nodes && row < nodes && col >= 0 && row >= 0){
210 if (row > col){i=row;row=col;col=i;}
212 else { E[k].v1=row; E[k].v2=col; k++;}
213 }
else if (gotone && k>=edges) {
214 printf(
"To many edges in file.\nLine %d, %s\n",line,thisline);
216 }
else if (gotone&&(col >= nodes || row >= nodes || col < 0 || row < 0)){
217 printf(
"Invalid node number.\nLine %d, %s\n",line,thisline);
222 *nnodes=nodes; *nedges=k;
229static int EdgeMatDestroy(
void*);
230static int EdgeMatView(
void*);
231static int EdgeMatVecVec(
void*,
double[],
int,
double *);
232static int EdgeMatDot(
void*,
double[],
int,
int,
double *);
233static int EdgeMatGetRank(
void*,
int*,
int);
234static int EdgeMatFactor(
void*);
235static int EdgeMatGetEig(
void*,
int,
double*,
double[],
int,
int[],
int*);
236static int EdgeMatAddRowMultiple(
void*,
int,
double,
double[],
int);
237static int EdgeMatAddMultiple(
void*,
double,
double[],
int,
int);
238static int EdgeMatGetRowNnz(
void*,
int,
int[],
int*,
int);
240static int EdgeMatDestroy(
void* AA){
243static int EdgeMatVecVec(
void* A,
double x[],
int n,
double *v){
244 EdgeMat*E=(EdgeMat*)A;
245 *v=2*x[E->v1]*x[E->v2];
248static int EdgeMatDot(
void* A,
double x[],
int nn,
int n,
double *v){
249 EdgeMat*E=(EdgeMat*)A;
250 int k=E->v2*(E->v2+1)/2 + E->v1;
254static int EdgeMatView(
void* A){
255 EdgeMat*E=(EdgeMat*)A;
256 printf(
" Row: %d, Column: %d\n",E->v1,E->v2);
257 printf(
" Row: %d, Column: %d\n",E->v2,E->v1);
260static int EdgeMatFactor(
void* A){
263static int EdgeMatGetRank(
void *A,
int*rank,
int n){
267static int EdgeMatGetEig(
void*A,
int neig,
double *eig,
double v[],
int n,
int spind[],
int *nind){
268 EdgeMat*E=(EdgeMat*)A;
269 double tt=1.0/sqrt(2.0);
270 memset((
void*)v,0,(n)*
sizeof(
double));
271 memset((
void*)spind,0,(n)*
sizeof(
int));
273 v[E->v1]=tt;v[E->v2]=tt;*eig=1;
274 spind[0]=E->v1;spind[1]=E->v2; *nind=2;
276 v[E->v1]=-tt;v[E->v2]=tt;*eig=-1;
277 spind[0]=E->v1;spind[1]=E->v2; *nind=2;
278 }
else { *eig=0;*nind=0;}
281static int EdgeMatGetRowNnz(
void*A,
int nrow,
int nz[],
int *nnzz,
int n){
282 EdgeMat*E=(EdgeMat*)A;
283 if (nrow==E->v1){ nz[E->v2]++; *nnzz=1;}
284 else if (nrow==E->v2){nz[E->v1]++; *nnzz=1;}
288static int EdgeMatAddRowMultiple(
void*A,
int nrow,
double dd,
double rrow[],
int n){
289 EdgeMat*E=(EdgeMat*)A;
290 if (nrow==E->v1){ rrow[E->v2]+=dd;}
291 else if (nrow==E->v2){rrow[E->v1]+=dd;}
294static int EdgeMatAddMultiple(
void*A,
double dd,
double vv[],
int nn,
int n){
295 EdgeMat*E=(EdgeMat*)A;
296 int k=E->v2*(E->v2+1)/2 + E->v1;
300static int EdgeMatFNorm(
void*A,
int n,
double *fnorm){
304static int EdgeMatCountNonzeros(
void*A,
int *nnz,
int n){
308static const char *datamatname=
"THETA EDGE MATRIX";
309static int EdgeMatOperationsInitialize(
struct DSDPDataMat_Ops* edgematoperator){
311 if (edgematoperator==NULL)
return 0;
313 edgematoperator->matfactor1=EdgeMatFactor;
314 edgematoperator->matgetrank=EdgeMatGetRank;
315 edgematoperator->matgeteig=EdgeMatGetEig;
316 edgematoperator->matvecvec=EdgeMatVecVec;
317 edgematoperator->matrownz=EdgeMatGetRowNnz;
318 edgematoperator->matdot=EdgeMatDot;
319 edgematoperator->matfnorm2=EdgeMatFNorm;
320 edgematoperator->matnnz=EdgeMatCountNonzeros;
321 edgematoperator->mataddrowmultiple=EdgeMatAddRowMultiple;
322 edgematoperator->mataddallmultiple=EdgeMatAddMultiple;
323 edgematoperator->matdestroy=EdgeMatDestroy;
324 edgematoperator->matview=EdgeMatView;
325 edgematoperator->matname=datamatname;
326 edgematoperator->id=25;
331static int OneMatDestroy(
void*);
332static int OneMatView(
void*);
333static int OneMatVecVec(
void*,
double[],
int,
double *);
334static int OneMatDot(
void*,
double[],
int,
int,
double *);
335static int OneMatGetRank(
void*,
int*,
int);
336static int OneMatFactor(
void*);
337static int OneMatGetEig(
void*,
int,
double*,
double[],
int,
int[],
int*);
338static int OneMatRowNnz(
void*,
int,
int[],
int*,
int);
339static int OneMatAddRowMultiple(
void*,
int,
double,
double[],
int);
340static int OneMatAddMultiple(
void*,
double,
double[],
int,
int);
343static int OneMatFactor(
void*A){
return 0;}
344static int OneMatGetRank(
void *A,
int *rank,
int n){*rank=1;
return 0;}
345static int OneMatFNorm2(
void*AA,
int n,
double *fnorm2){*fnorm2=1.0*n*n;
return 0;}
346static int OneMatCountNonzeros(
void*AA,
int *nnz,
int n){*nnz=n*n;
return 0;}
347static int OneMatDot(
void* A,
double x[],
int nn,
int n,
double *v){
351 for (j=0;j<i;j++,x++){dtmp+= (*x);}
357static int OneMatVecVec(
void* A,
double x[],
int n,
double *v){
360 for (i=0; i<n; i++){dtmp+=x[i];}
364static int OneMatAddMultiple(
void*A,
double ddd,
double vv[],
int nn,
int n){
367 for (j=0;j<i;j++,vv++){(*vv)+=-ddd;}
372static int OneMatAddRowMultiple(
void*A,
int nrow,
double ddd,
double row[],
int n){
374 for (i=0;i<n;i++){row[i] -= ddd;}
378static int OneMatGetEig(
void*A,
int neig,
double *eig,
double v[],
int n,
int spind[],
int *nind){
380 if (neig==0){ *eig=-1;
for (i=0;i<n;i++){ v[i]=1.0; spind[i]=i;} *nind=n;
381 }
else { *eig=0;
for (i=0;i<n;i++){ v[i]=0.0; } *nind=0;
385static int OneMatRowNnz(
void*A,
int row,
int nz[],
int *nnz,
int n){
387 for (i=0;i<n;i++){ nz[i]++; }
391static int OneMatView(
void* AA){
392 printf(
"Every element of the matrix is the same: -1\n");
395static int OneMatDestroy(
void* A){
return 0;}
397static const char *mat1name=
"THETA ALL ELEMENTS EQUAL -ONE";
400 if (mat1ops==NULL)
return 0;
402 mat1ops->matfactor1=OneMatFactor;
403 mat1ops->matgetrank=OneMatGetRank;
404 mat1ops->matgeteig=OneMatGetEig;
405 mat1ops->matvecvec=OneMatVecVec;
406 mat1ops->matdot=OneMatDot;
407 mat1ops->mataddrowmultiple=OneMatAddRowMultiple;
408 mat1ops->mataddallmultiple=OneMatAddMultiple;
409 mat1ops->matdestroy=OneMatDestroy;
410 mat1ops->matview=OneMatView;
411 mat1ops->matrownz=OneMatRowNnz;
412 mat1ops->matfnorm2=OneMatFNorm2;
413 mat1ops->matnnz=OneMatCountNonzeros;
415 mat1ops->matname=mat1name;
420static int IdentityMatDestroy(
void*);
421static int IdentityMatView(
void*);
422static int IdentityMatVecVec(
void*,
double[],
int,
double *);
423static int IdentityMatDot(
void*,
double[],
int,
int,
double *);
424static int IdentityMatGetRank(
void*,
int*,
int);
425static int IdentityMatFactor(
void*);
426static int IdentityMatGetEig(
void*,
int,
double*,
double[],
int,
int[],
int*);
427static int IdentityMatAddRowMultiple(
void*,
int,
double,
double[],
int);
428static int IdentityMatAddMultiple(
void*,
double,
double[],
int,
int);
429static int IdentityMatGetRowNnz(
void*,
int,
int[],
int*,
int);
431static int IdentityMatDestroy(
void* AA){
return 0;}
432static int IdentityMatFNorm2(
void* AA,
int n,
double *v){*v=1.0*n;
return 0;}
433static int IdentityMatGetRank(
void *AA,
int*rank,
int n){*rank=n;
return 0;}
434static int IdentityMatFactor(
void*A){
return 0;}
435static int IdentityMatCountNonzeros(
void*A,
int *nnz,
int n){*nnz=n;
return 0;}
436static int IdentityMatVecVec(
void* AA,
double x[],
int n,
double *v){
439 for (i=0;i<n;i++){ *v+=x[i]*x[i]; }
442static int IdentityMatDot(
void* AA,
double x[],
int nn,
int n,
double *v){
445 for (i=0;i<n;i++){ vv+=x[((i+1)*(i+2))/2-1];}
449static int IdentityMatView(
void* AA){
450 printf(
"Identity matrix: All Diagonal elements equal 1.0\n");
453static int IdentityMatGetEig(
void*AA,
int neig,
double *eig,
double v[],
int n,
int spind[],
int *nind){
454 if (neig<0 || neig>=n){ *eig=0;
return 0;}
455 memset((
void*)v,0,(n)*
sizeof(
double));
456 *eig=1.0; v[neig]=1.0; spind[0]=neig; *nind=1;
459static int IdentityMatGetRowNnz(
void*A,
int nrow,
int nz[],
int *nnzz,
int n){
460 if (nrow>=0 && nrow < n){
466static int IdentityMatAddRowMultiple(
void*A,
int nrow,
double dd,
double rrow[],
int n){
467 rrow[nrow] += dd;
return 0;
469static int IdentityMatAddMultiple(
void*A,
double dd,
double vv[],
int nn,
int n){
471 for (i=0;i<n;i++){ vv[(i+1)*(i+2)/2-1] += dd;}
475static const char *eyematname=
"THETA IDENTITY MATRIX";
476static int IdentitymatOperationsInitialize(
struct DSDPDataMat_Ops* imatops){
478 if (imatops==NULL)
return 0;
480 imatops->matfactor1=IdentityMatFactor;
481 imatops->matgetrank=IdentityMatGetRank;
482 imatops->matgeteig=IdentityMatGetEig;
483 imatops->matvecvec=IdentityMatVecVec;
484 imatops->matrownz=IdentityMatGetRowNnz;
485 imatops->matdot=IdentityMatDot;
486 imatops->matfnorm2=IdentityMatFNorm2;
487 imatops->matnnz=IdentityMatCountNonzeros;
488 imatops->mataddrowmultiple=IdentityMatAddRowMultiple;
489 imatops->mataddallmultiple=IdentityMatAddMultiple;
490 imatops->matdestroy=IdentityMatDestroy;
491 imatops->matview=IdentityMatView;
493 imatops->matname=eyematname;
The API to DSDP for those applications using DSDP as a subroutine library.
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Table of function pointers that operate on the data matrix.
Internal structures for the DSDP solver.
Internal structure for semidefinite cone.