1 #include "dsdp5.h"
2 #include <string.h>
3 #include <math.h>
4 #include <stdio.h>
5 #include <stdlib.h>
6
7 /*! \file color.c
8 \brief Second Basic Example: Read graph from file, formulate the SDP relaxation of k-coloring
9 problem, solve using DSDP, and apply randomized algorithm to generate
10 approximate solutions.
11 */
12 char help2[]="\nA positive semidefinite relaxation of the\n\
13 graph k coloring problem can be rewritten as\n\n\
14 Find X>=0 \n\
15 such that X_ij <= 1 - 1/(k-1) for all edges (i,j).\n\
16 ";
17
18 char help[]="DSDP Usage: color <graph filename> ";
19
20 static int ReadGraph(char*,int *, int *,int**, int **, double **);
21 static int ParseGraphline(char*,int*,int*,double*,int*);
22 static int RandomizedColor(DSDP, SDPCone, int, int[], int[], int);
23 int MinColoring(int argc,char *argv[]);
24
25
main(int argc,char * argv[])26 int main(int argc,char *argv[]){
27 int info;
28 info=MinColoring(argc,argv);
29 return 0;
30 }
31
32 /*!
33 \fn int MinColoring(int argc,char *argv[]);
34 \param argc number of command line arguments
35 \param argv command line arguments
36 \ingroup Examples
37 \brief SDP relaxation of k-coloring problem
38 */
MinColoring(int argc,char * argv[])39 int MinColoring(int argc,char *argv[]){
40
41 int i,kk,vari,info;
42 int *node1,*node2,nedges,nnodes;
43 int *iptr1,*iptr2;
44 double *weight,*yy1,*yy2,bb;
45 DSDPTerminationReason reason;
46 SDPCone sdpcone;
47 BCone bcone;
48 DSDP dsdp;
49
50 if (argc<2){ printf("%s\n%s",help2,help); return(1); }
51
52 info = ReadGraph(argv[1],&nnodes,&nedges,&node1,&node2,&weight);
53 if (info){ printf("Problem reading file\n"); return 1; }
54
55 info = DSDPCreate(nnodes+nedges,&dsdp);
56
57 info = DSDPCreateSDPCone(dsdp,1,&sdpcone);
58 info = SDPConeSetBlockSize(sdpcone,0,nnodes);
59 info = SDPConeSetSparsity(sdpcone,0,nnodes+nedges+1);
60
61 info = DSDPCreateBCone(dsdp,&bcone);
62
63 if (info){ printf("Out of memory\n"); return 1; }
64
65
66 /* Formulate the problem from the data */
67 /* Create data matrices */
68
69 /* Diagonal elements of X(i,i) must equal 1.0 */
70 iptr1=(int*)malloc(nnodes*sizeof(int));
71 yy1=(double*)malloc(nnodes*sizeof(double));
72 for (i=0;i<nnodes;i++){
73 iptr1[i]=(i+2)*(i+1)/2-1;
74 yy1[i]=1.0;
75 }
76 for (i=0;i<nnodes;i++){
77 info=SDPConeSetASparseVecMat(sdpcone,0,i+1,nnodes,1.0,0,iptr1+i,yy1+i,1);
78 if (info) printf("ERROR 1: %d \n",i);
79 info=DSDPSetDualObjective(dsdp,i+1,1.0);
80 }
81
82 /* For each nonzero element (i,j) of the matrix, X(i,j) must be less than 1 - 1/nnodes */
83 bb=2-2.0/nnodes;
84 iptr2=(int*)malloc(nedges*sizeof(int));
85 yy2=(double*)malloc(nedges*sizeof(double));
86 for (i=0;i<nedges;i++){
87 iptr2[i]=(node1[i])*(node1[i]+1)/2+node2[i];
88 yy2[i]=1.0;
89 }
90 info = BConeAllocateBounds(bcone,nedges);
91 for (i=0; i<nedges; i++){
92 vari=nnodes+i+1;
93 info = SDPConeSetSparseVecMat(sdpcone,0,vari,nnodes,0,iptr2+i,yy2+i,1);
94 if (info) printf("ERROR 2: %d %d \n",i,vari);
95 info = BConeSetPSlackVariable(bcone,vari);
96 if (info) printf("ERROR 3: %d %d \n",i,vari);
97 info = DSDPSetDualObjective(dsdp,vari,bb);
98 }
99
100
101 /* Get read to go */
102 info=DSDPSetPotentialParameter(dsdp,5);
103
104 for (kk=1; kk<argc-1; kk++){
105 if (strncmp(argv[kk],"-dloginfo",8)==0){
106 info=DSDPLogInfoAllow(DSDP_TRUE,0);
107 } else if (strncmp(argv[kk],"-params",7)==0){
108 info=DSDPReadOptions(dsdp,argv[kk+1]);
109 } else if (strncmp(argv[kk],"-help",7)==0){
110 printf("%s\n",help);
111 }
112 }
113 info=DSDPSetOptions(dsdp,argv,argc);
114
115 if (info){ printf("Out of memory\n"); return 1; }
116 info = DSDPSetStandardMonitor(dsdp,1);
117
118 info = DSDPSetup(dsdp);
119 if (info){ printf("Out of memory\n"); return 1; }
120
121 info = DSDPSolve(dsdp);
122 if (info){ printf("Numerical error\n"); return 1; }
123 info = DSDPStopReason(dsdp,&reason);
124
125 if (reason!=DSDP_INFEASIBLE_START){ /* Randomized solution strategy */
126 info=RandomizedColor(dsdp, sdpcone, nnodes, node1, node2, nedges);
127 }
128
129 info = DSDPDestroy(dsdp);
130
131 free(node1);free(node2);free(weight);
132 free(iptr1);
133 free(yy1);
134 free(iptr2);
135 free(yy2);
136
137 return 0;
138 }
139
GetXRow(double xmat[],double xrow[],int row,int n)140 static int GetXRow(double xmat[],double xrow[],int row,int n){
141 int i,i1=row*(row+1)/2;
142 for (i=0;i<row;i++){xrow[i]=xmat[i1+i];}
143 for (i=row;i<n;i++){xrow[i]=xmat[i1+row];i1+=i+1;}
144 return 0;
145 }
146
147 typedef struct {
148 int index;double val;
149 } orderVec;
150
cut_comp(const void * e1,const void * e2)151 static int cut_comp(const void *e1,const void *e2){ /* Used in qsort routine */
152 double d1=((orderVec*)e1)->val, d2=((orderVec*)e2)->val;
153 if (d1<d2) return (1);
154 else if (d1>d2) return (-1);
155 return(0);
156 }
157
RemoveNode(int node,int node1[],int node2[],int * nedges)158 static int RemoveNode(int node, int node1[], int node2[], int *nedges){
159 int i,nnedges=*nedges;
160 for (i=0;i<nnedges;i++){
161 if (node1[i]==node || node2[i]==node){
162 node1[i]=node1[nnedges-1];
163 node2[i]=node2[nnedges-1];
164 nnedges--;
165 if (i < nnedges) i--;
166 }
167 }
168 *nedges=nnedges;
169 return 0;
170 }
171
Connected(int n1,int n2,int node1[],int node2[],int nedges)172 static int Connected(int n1, int n2, int node1[], int node2[], int nedges){
173 int i;
174 if (n1==n2) return 1;
175 for (i=0;i<nedges;i++){
176 if (node1[i]==n1 && node2[i]==n2){ return 1;}
177 if (node1[i]==n2 && node2[i]==n1){ return 1;}
178 }
179 return 0;
180 }
181
HighDegree(int node1[],int node2[],int nedges,int iwork[],int nnodes)182 static int HighDegree(int node1[], int node2[], int nedges, int iwork[], int nnodes){
183 int i,nmax=0,maxdegree=-1;
184 for (i=0;i<nnodes;i++) iwork[i]=0;
185 for (i=0;i<nedges;i++){
186 iwork[node1[i]]++; iwork[node2[i]]++;
187 }
188 for (i=0;i<nnodes;i++){ if (iwork[i]>maxdegree){nmax=i; maxdegree=iwork[i];} }
189 return nmax;
190 }
191
First(int coloring[],int nnodes)192 static int First(int coloring[], int nnodes){
193 int i,nmax=nnodes;
194 for (i=0;i<nnodes;i++){
195 if (coloring[i]==0){
196 nmax=i; return nmax;
197 }
198 }
199 return -1;
200 }
201
RandomizedColor(DSDP dsdp,SDPCone sdpcone,int nnodes,int node1[],int node2[],int nedges)202 static int RandomizedColor(DSDP dsdp, SDPCone sdpcone, int nnodes, int node1[], int node2[], int nedges){
203 int i,j,nodek,nn,info,flag,coloring=0,maxvertex;
204 int *degree,*color,*cgroup,ngsize,uncolored=nnodes;
205 int tnedges=nedges;
206 double *xrow,*xptr;
207 orderVec *vorder;
208
209 xrow=(double*)malloc(nnodes*sizeof(double));
210 color=(int*)malloc(nnodes*sizeof(int));
211 cgroup=(int*)malloc(nnodes*sizeof(int));
212 degree=(int*)malloc(nnodes*sizeof(int));
213 vorder=(orderVec*)malloc(nnodes*sizeof(orderVec));
214
215 for (i=0;i<nnodes;i++){ color[i]=0;}
216 info=DSDPComputeX(dsdp);
217 info=SDPConeGetXArray(sdpcone,0,&xptr,&nn);
218
219 while (uncolored>0){
220
221 coloring++;
222
223 maxvertex=First(color,nnodes);
224 maxvertex=HighDegree(node1,node2,tnedges,degree,nnodes);
225
226 cgroup[0]=maxvertex;ngsize=1;
227
228 info=GetXRow(xptr,xrow,maxvertex,nnodes);
229
230 for (i=0;i<nnodes;i++){vorder[i].index=i; vorder[i].val = xrow[i];}
231 qsort( (void*)vorder, nnodes, sizeof(orderVec), cut_comp);
232
233 for (i=0;i<nnodes;i++){
234 nodek=vorder[i].index;
235 if (color[nodek]==0){
236 for (flag=0,j=0;j<ngsize;j++){
237 if (Connected(nodek,cgroup[j],node1,node2,tnedges) ){flag=1;}
238 }
239 if (flag==0){ cgroup[ngsize]=nodek; ngsize++; }
240 }
241 }
242 for (i=0;i<ngsize;i++){
243 color[cgroup[i]]=coloring; uncolored--;
244 RemoveNode(cgroup[i],node1,node2,&tnedges);
245 }
246 }
247 printf("\nCOLORS: %d\n",coloring);
248 free(xrow);
249 free(color);
250 free(cgroup);
251 free(degree);
252 free(vorder);
253 return(0);
254 }
255
256
257 #define BUFFERSIZ 100
258
259 #undef __FUNCT__
260 #define __FUNCT__ "ParseGraphline"
ParseGraphline(char * thisline,int * row,int * col,double * value,int * gotem)261 static int ParseGraphline(char * thisline,int *row,int *col,double *value,
262 int *gotem){
263
264 int temp;
265 int rtmp,coltmp;
266
267 rtmp=-1, coltmp=-1, *value=0.0;
268 temp=sscanf(thisline,"%d %d %lf",&rtmp,&coltmp,value);
269 if (temp==3 && rtmp>0 && coltmp>0) *gotem=3;
270 else if (temp==2 && rtmp>0 && coltmp>0){ *value = 1.0; *gotem=3;}
271 else *gotem=0;
272 *row=rtmp-1; *col=coltmp-1;
273
274 return(0);
275 }
276
277
278 #undef __FUNCT__
279 #define __FUNCT__ "ReadGraph"
ReadGraph(char * filename,int * nnodes,int * nedges,int ** n1,int ** n2,double ** wght)280 int ReadGraph(char* filename,int *nnodes, int *nedges,
281 int**n1, int ** n2, double **wght){
282
283 FILE*fp;
284 char thisline[BUFFERSIZ]="*";
285 int i,k=0,line=0,nodes,edges,gotone=3;
286 int *node1,*node2;
287 double *weight;
288 int info,row,col;
289 double value;
290
291 fp=fopen(filename,"r");
292 if (!fp){printf("Cannot open file %s !",filename);return(1);}
293
294 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){
295 fgets(thisline,BUFFERSIZ,fp); line++;
296 }
297
298 if (sscanf(thisline,"%d %d",&nodes, &edges)!=2){
299 printf("First line must contain the number of nodes and number of edges\n");
300 return 1;
301 }
302
303 node1=(int*)malloc(edges*sizeof(int));
304 node2=(int*)malloc(edges*sizeof(int));
305 weight=(double*)malloc(edges*sizeof(double));
306
307 for (i=0; i<edges; i++){ node1[i]=0;node2[i]=0;weight[i]=0.0;}
308
309 while(!feof(fp) && gotone){
310 thisline[0]='\0';
311 fgets(thisline,BUFFERSIZ,fp); line++;
312 info = ParseGraphline(thisline,&row,&col,&value,&gotone);
313 if (gotone && value!=0.0 && k<edges &&
314 col < nodes && row < nodes && col >= 0 && row >= 0){
315 if (row<col){info=row;row=col;col=info;}
316 if (row == col){}
317 else {
318 node1[k]=row; node2[k]=col;
319 weight[k]=value; k++;
320 }
321 } else if (gotone && k>=edges) {
322 printf("To many edges in file.\nLine %d, %s\n",line,thisline);
323 return 1;
324 } else if (gotone&&(col >= nodes || row >= nodes || col < 0 || row < 0)){
325 printf("Invalid node number.\nLine %d, %s\n",line,thisline);
326 return 1;
327 }
328 }
329 *nnodes=nodes; *nedges=edges;
330 *n1=node1; *n2=node2; *wght=weight;
331 return 0;
332 }
333