DSDP
theta.c
Go to the documentation of this file.
1#include "dsdp5.h"
2#include <stdio.h>
3#include <string.h>
4#include <math.h>
5#include <stdlib.h>
6
11char help[]="\n\
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";
19
20typedef struct {
21 int v1,v2;
22} EdgeMat;
23
24#include "../src/sdp/dsdpdatamat_impl.h"
25extern int SDPConeAddDataMatrix(SDPCone,int, int, int, char, struct DSDPDataMat_Ops*, void*);
26
27
28/* These variables and routines are for customized SDP Data Matrices */
29static struct DSDPDataMat_Ops edgematop;
30static int EdgeMatOperationsInitialize(struct DSDPDataMat_Ops*);
31
32static struct DSDPDataMat_Ops onematops;
33static int OneMatOpsInitialize(struct DSDPDataMat_Ops*);
34
35static struct DSDPDataMat_Ops identitymatops;
36static int IdentitymatOperationsInitialize(struct DSDPDataMat_Ops*);
37
38static int ReadGraphFromFile(char*,int*, int*, EdgeMat*[]);
39int SetThetaData(DSDP, SDPCone, int, int, EdgeMat[]);
40int LovaszTheta(int argc,char *argv[]);
41
42int main(int argc,char *argv[]){
43 int info;
44 info=LovaszTheta(argc,argv);
45 return 0;
46}
47
56int LovaszTheta(int argc,char *argv[]){
57
58 int info,kk,nedges,nnodes;
59 EdgeMat *Edges;
60 DSDP dsdp;
61 SDPCone sdpcone;
62
63 if (argc<2){ printf("%s",help); return(1); }
64
65 info = ReadGraphFromFile(argv[1],&nnodes,&nedges,&Edges);
66 if (info){ printf("Problem reading file\n"); return 1; }
67
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; }
76
77 info = DSDPSetGapTolerance(dsdp,0.001);
78 info = DSDPSetZBar(dsdp,nnodes+1);
79 info = DSDPReuseMatrix(dsdp,10);
80
81 for (kk=1; kk<argc-1; kk++){
82 if (strncmp(argv[kk],"-dloginfo",8)==0){
83 info=DSDPLogInfoAllow(DSDP_TRUE,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){
87 printf("%s\n",help);
88 }
89 }
90 info=DSDPSetOptions(dsdp,argv,argc);
91
92 if (info){ printf("Out of memory\n"); return 1; }
93 info = DSDPSetStandardMonitor(dsdp,1);
94 if (info){ printf("Monitor Problem \n"); return 1; }
95
96 info = DSDPSetup(dsdp);
97 if (info){ printf("Out of memory\n"); return 1; }
98 if (0==1){info=SDPConeCheckData(sdpcone);}
99
100 info = DSDPSolve(dsdp);
101 if (info){ printf("Numberical error\n"); return 1; }
102 info=DSDPComputeX(dsdp);DSDPCHKERR(info);
103
104 if (0==1){ /* Look at the solution */
105 int n; double *xx;
106 info=SDPConeGetXArray(sdpcone,0,&xx,&n);
107 }
108
109 info = DSDPDestroy(dsdp);
110 free(Edges);
111
112 return 0;
113} /* main */
114
126int SetThetaData(DSDP dsdp, SDPCone sdpcone, int nodes, int edges, EdgeMat Edge[]){
127
128 int i,info;
129
130 /* Create data matrices - these are all custom types */
131
132 /* The c matrix has all elements equal to 1.0 */
133 info=OneMatOpsInitialize(&onematops);
134 info=SDPConeAddDataMatrix(sdpcone,0,0,nodes,'P',&onematops,0);
135
136 /* For each edge connecting nodes i and j, X(i,j)=X(j,i)=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;
143 }
144
145 /* The trace of X must equal 1.0 */
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);
150
151 /* The initial point y is feasible and near the central path */
152 info = DSDPSetR0(dsdp,0);
153
154 return(0);
155}
156
157#define BUFFERSIZ 100
158
159#undef __FUNCT__
160#define __FUNCT__ "ParseGraphline"
161static int ParseGraphline(char * thisline,int *row,int *col,double *value,
162 int *gotem){
163
164 int temp;
165 int rtmp,coltmp;
166
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;}
171 else *gotem=0;
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);
175 }
176 return(0);
177}
178
179#undef __FUNCT__
180#define __FUNCT__ "ReadGraphFromFile"
181static int ReadGraphFromFile(char* filename,int *nnodes, int *nedges, EdgeMat**EE){
182
183 FILE*fp;
184 char thisline[BUFFERSIZ]="*";
185 int i,k=0,line=0,nodes,edges,gotone=3;
186 int info,row,col;
187 double value;
188 EdgeMat *E;
189
190 fp=fopen(filename,"r");
191 if (!fp){ printf("Cannot open file %s !",filename); return(1); }
192
193 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){
194 fgets(thisline,BUFFERSIZ,fp); line++; }
195
196 if (sscanf(thisline,"%d %d",&nodes, &edges)!=2){
197 printf("First line must contain the number of nodes and number of edges\n");
198 return 1;
199 }
200
201 E=(EdgeMat*)malloc(edges*sizeof(EdgeMat)); *EE=E;
202 for (i=0; i<edges; i++){ E[i].v1=0; E[i].v2=0; }
203
204 while(!feof(fp) && gotone){
205 thisline[0]='\0';
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;}
211 if (row == col){}
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);
215 return 1;
216 } else if (gotone&&(col >= nodes || row >= nodes || col < 0 || row < 0)){
217 printf("Invalid node number.\nLine %d, %s\n",line,thisline);
218 return 1;
219 }
220 }
221
222 *nnodes=nodes; *nedges=k;
223
224 return 0;
225}
226
227
228/* SPecial Matrix where each edge represents a constraint matrix */
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);
239
240static int EdgeMatDestroy(void* AA){
241 return 0;
242}
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];
246 return 0;
247}
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;
251 *v=2*x[k];
252 return 0;
253}
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);
258 return 0;
259}
260static int EdgeMatFactor(void* A){
261 return 0;
262}
263static int EdgeMatGetRank(void *A, int*rank, int n){
264 *rank=2;
265 return 0;
266}
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));
272 if (neig==0){
273 v[E->v1]=tt;v[E->v2]=tt;*eig=1;
274 spind[0]=E->v1;spind[1]=E->v2; *nind=2;
275 } else if (neig==1){
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;}
279 return 0;
280}
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;}
285 else {*nnzz=0;}
286 return 0;
287}
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;}
292 return 0;
293}
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;
297 vv[k]+=dd;
298 return 0;
299}
300static int EdgeMatFNorm(void*A, int n, double *fnorm){
301 *fnorm=2.0;
302 return 0;
303}
304static int EdgeMatCountNonzeros(void*A, int *nnz, int n){
305 *nnz=1;
306 return 0;
307}
308static const char *datamatname="THETA EDGE MATRIX";
309static int EdgeMatOperationsInitialize(struct DSDPDataMat_Ops* edgematoperator){
310 int info;
311 if (edgematoperator==NULL) return 0;
312 info=DSDPDataMatOpsInitialize(edgematoperator); if (info){ return info;}
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;
327 return 0;
328}
329
330/* SPecial Matrix where all elements equal negative one */
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);
341
342
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){
348 double dtmp=0.0;
349 int i,j;
350 for (i=0;i<n;i++){
351 for (j=0;j<i;j++,x++){dtmp+= (*x);}
352 dtmp+= (*x);x++;
353 }
354 *v=-2*dtmp;
355 return 0;
356}
357static int OneMatVecVec(void* A, double x[], int n, double *v){
358 double dtmp=0.0;
359 int i;
360 for (i=0; i<n; i++){dtmp+=x[i];}
361 *v=-dtmp*dtmp;
362 return 0;
363}
364static int OneMatAddMultiple(void*A, double ddd, double vv[], int nn, int n){
365 int i,j;
366 for (i=0;i<n;i++){
367 for (j=0;j<i;j++,vv++){(*vv)+=-ddd;}
368 (*vv)+=-ddd; vv++;
369 }
370 return 0;
371}
372static int OneMatAddRowMultiple(void*A, int nrow, double ddd, double row[], int n){
373 int i;
374 for (i=0;i<n;i++){row[i] -= ddd;}
375 row[nrow] -= ddd;
376 return 0;
377}
378static int OneMatGetEig(void*A, int neig, double *eig, double v[], int n, int spind[], int *nind){
379 int i;
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;
382 }
383 return 0;
384}
385static int OneMatRowNnz(void*A, int row, int nz[], int *nnz, int n){
386 int i;
387 for (i=0;i<n;i++){ nz[i]++; }
388 *nnz=n;
389 return 0;
390}
391static int OneMatView(void* AA){
392 printf("Every element of the matrix is the same: -1\n");
393 return 0;
394}
395static int OneMatDestroy(void* A){return 0;}
396
397static const char *mat1name="THETA ALL ELEMENTS EQUAL -ONE";
398static int OneMatOpsInitialize(struct DSDPDataMat_Ops* mat1ops){
399 int info;
400 if (mat1ops==NULL) return 0;
401 info=DSDPDataMatOpsInitialize(mat1ops); DSDPCHKERR(info);
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;
414 mat1ops->id=18;
415 mat1ops->matname=mat1name;
416 return 0;
417}
418
419/* Special Matrix for the Identity */
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);
430
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){
437 int i;
438 *v=0;
439 for (i=0;i<n;i++){ *v+=x[i]*x[i]; }
440 return 0;
441}
442static int IdentityMatDot(void* AA, double x[], int nn, int n, double *v){
443 int i;
444 double vv=0;
445 for (i=0;i<n;i++){ vv+=x[((i+1)*(i+2))/2-1];}
446 *v = 2*vv;
447 return 0;
448}
449static int IdentityMatView(void* AA){
450 printf("Identity matrix: All Diagonal elements equal 1.0\n");
451 return 0;
452}
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;
457 return 0;
458}
459static int IdentityMatGetRowNnz(void*A, int nrow, int nz[], int *nnzz, int n){
460 if (nrow>=0 && nrow < n){
461 *nnzz=1; nz[nrow]++;
462 } else { *nnzz=0;
463 }
464 return 0;
465}
466static int IdentityMatAddRowMultiple(void*A, int nrow, double dd, double rrow[], int n){
467 rrow[nrow] += dd;return 0;
468}
469static int IdentityMatAddMultiple(void*A, double dd, double vv[], int nn, int n){
470 int i;
471 for (i=0;i<n;i++){ vv[(i+1)*(i+2)/2-1] += dd;}
472 return 0;
473}
474
475static const char *eyematname="THETA IDENTITY MATRIX";
476static int IdentitymatOperationsInitialize(struct DSDPDataMat_Ops* imatops){
477 int info;
478 if (imatops==NULL) return 0;
479 info=DSDPDataMatOpsInitialize(imatops); if (info){return info;}
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;
492 imatops->id=12;
493 imatops->matname=eyematname;
494 return 0;
495}
The API to DSDP for those applications using DSDP as a subroutine library.
@ DSDP_TRUE
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Definition dsdpdatamat.c:47
Table of function pointers that operate on the data matrix.
Internal structures for the DSDP solver.
Definition dsdp.h:65
Internal structure for semidefinite cone.
Definition dsdpsdp.h:80