1 #include "dsdp5.h"
2 #include <string.h>
3 #include <math.h>
4 #include <stdio.h>
5 #include <stdlib.h>
6 /*! \file maxcut.c
7 \brief Most Basic Example: read graph from file, formulate the SDP relaxation of maximum cut
8 problem, solve using DSDP, and apply randomized algorithm to generate
9 approximate solutions.
10 */
11
12 char help[]="\
13 DSDP Usage: maxcut <graph filename> \n\
14 -gaptol <relative duality gap: default is 0.001> \n\
15 -maxit <maximum iterates: default is 200> \n";
16
17 static int ReadGraph(char*,int *, int *,int**, int **, double **);
18 static int TCheckArgs(DSDP,int,char **);
19 static int ParseGraphline(char*,int*,int*,double*,int*);
20 int MaxCutRandomized(SDPCone,int);
21 int MaxCut(int,int, int[], int[], double[]);
22
23
main(int argc,char * argv[])24 int main(int argc,char *argv[]){
25 int info;
26 int *node1,*node2,nedges,nnodes;
27 double *weight;
28
29 if (argc<2){ printf("%s",help); return(1); }
30
31 info = ReadGraph(argv[1],&nnodes,&nedges,&node1,&node2,&weight);
32 if (info){ printf("Problem reading file\n"); return 1; }
33
34 MaxCut(nnodes,nedges,node1,node2,weight);
35
36 free(node1);free(node2);free(weight);
37 return 0;
38 }
39
40 /*!
41 \fn int MaxCut(int nnodes,int nedged, int node1[], int node2[], double weight[]);
42 \brief Formulate and solve the SDP relaxation of the Maximum Cut problem
43 \ingroup Examples
44 \param nnodes number of nodes in graph
45 \param nedges number of edges in graph
46 \param node1 first node of each edge
47 \param node2 second node of each edge
48 \param weight weight of each edge
49 \note This routine is an example! It is not part of the solver library.
50 */
MaxCut(int nnodes,int nedges,int node1[],int node2[],double weight[])51 int MaxCut(int nnodes,int nedges, int node1[], int node2[], double weight[]){
52
53 int i,info;
54 int *indd,*iptr;
55 double *yy,*val,*diag,tval=0;
56 DSDPTerminationReason reason;
57 SDPCone sdpcone;
58 DSDP dsdp;
59
60 info = DSDPCreate(nnodes,&dsdp);
61 info = DSDPCreateSDPCone(dsdp,1,&sdpcone);
62
63 if (info){ printf("Out of memory\n"); return 1; }
64
65 info = SDPConeSetBlockSize(sdpcone,0,nnodes);
66
67
68 /* Formulate the problem from the data */
69 /*
70 Diagonal elements equal 1.0
71 Create Constraint matrix A_i for i=1, ..., nnodes.
72 that has a single nonzero element.
73 */
74 diag=(double*)malloc(nnodes*sizeof(double));
75 iptr=(int*)malloc(nnodes*sizeof(int));
76 for (i=0;i<nnodes;i++){
77 iptr[i]=i*(i+1)/2+i;
78 diag[i]=1.0;
79 }
80
81 for (i=0;i<nnodes;i++){
82 info = DSDPSetDualObjective(dsdp,i+1,1.0);
83 info = SDPConeSetASparseVecMat(sdpcone,0,i+1,nnodes,1.0,0,iptr+i,diag+i,1);
84 if (0==1){
85 printf("Matrix: %d\n",i+1);
86 info = SDPConeViewDataMatrix(sdpcone,0,i+1);
87 }
88 }
89
90 /* C matrix is the Laplacian of the adjacency matrix */
91 /* Also compute a feasible initial point y such that S>=0 */
92 yy=(double*)malloc(nnodes*sizeof(double));
93 for (i=0;i<nnodes;i++){yy[i]=0.0;}
94 indd=(int*)malloc((nnodes+nedges)*sizeof(int));
95 val=(double*)malloc((nnodes+nedges)*sizeof(double));
96 for (i=0;i<nnodes+nedges;i++){indd[i]=0;}
97 for (i=0;i<nnodes;i++){indd[nedges+i]=i*(i+1)/2+i;}
98 for (i=0;i<nnodes+nedges;i++){val[i]=0;}
99 for (i=0;i<nedges;i++){
100 indd[i]=(node1[i])*(node1[i]+1)/2 + node2[i];
101 tval+=fabs(weight[i]);
102 val[i]=weight[i]/4;
103 val[nedges+node1[i]]-=weight[i]/4;
104 val[nedges+node2[i]]-=weight[i]/4;
105 yy[node1[i]]-= fabs(weight[i]/2);
106 yy[node2[i]]-= fabs(weight[i]/2);
107 }
108
109 if (0){
110 info = SDPConeSetASparseVecMat(sdpcone,0,0,nnodes,1.0,0,indd,val,nedges+nnodes);
111 } else { /* Equivalent */
112 info = SDPConeSetASparseVecMat(sdpcone,0,0,nnodes,1.0,0,indd,val,nedges);
113 info = SDPConeAddASparseVecMat(sdpcone,0,0,nnodes,1.0,0,indd+nedges,val+nedges,nnodes);
114 }
115 if (0==1){ info = SDPConeViewDataMatrix(sdpcone,0,0);}
116
117 /* Initial Point */
118 info = DSDPSetR0(dsdp,0.0);
119 info = DSDPSetZBar(dsdp,10*tval+1.0);
120 for (i=0;i<nnodes; i++){
121 info = DSDPSetY0(dsdp,i+1,10*yy[i]);
122 }
123 if (info) return info;
124
125 /* Get read to go */
126 info=DSDPSetGapTolerance(dsdp,0.001);
127 info=DSDPSetPotentialParameter(dsdp,5);
128 info=DSDPReuseMatrix(dsdp,0);
129 info=DSDPSetPNormTolerance(dsdp,1.0);
130 /*
131 info = TCheckArgs(dsdp,argc,argv);
132 */
133
134 if (info){ printf("Out of memory\n"); return 1; }
135 info = DSDPSetStandardMonitor(dsdp,1);
136
137 info = DSDPSetup(dsdp);
138 if (info){ printf("Out of memory\n"); return 1; }
139
140 info = DSDPSolve(dsdp);
141 if (info){ printf("Numerical error\n"); return 1; }
142 info = DSDPStopReason(dsdp,&reason);
143
144 if (reason!=DSDP_INFEASIBLE_START){ /* Randomized solution strategy */
145 info=MaxCutRandomized(sdpcone,nnodes);
146 if (0==1){ /* Look at the solution */
147 int n; double *xx,*y=diag;
148 info=DSDPGetY(dsdp,y,nnodes);
149 info=DSDPComputeX(dsdp);DSDPCHKERR(info);
150 info=SDPConeGetXArray(sdpcone,0,&xx,&n);
151 }
152 }
153 info = DSDPDestroy(dsdp);
154
155 free(iptr);
156 free(yy);
157 free(indd);
158 free(val);
159 free(diag);
160
161 return 0;
162 } /* main */
163
164
165
166 /*!
167 int MaxCutRandomized(SDPCone sdpcone,int nnodes);
168 \brief Apply the Goemens and Williamson randomized cut algorithm to the SDP relaxation of the max-cut problem
169 \param sdpcone the SDP cone
170 \param nnodes number of nodes in the graph
171 \note This routine is an example! It is not part of the solver library.
172 \ingroup Examples
173 \sa MaxCut()
174 */
MaxCutRandomized(SDPCone sdpcone,int nnodes)175 int MaxCutRandomized(SDPCone sdpcone,int nnodes){
176 int i,j,derror,info;
177 double dd,scal=RAND_MAX,*vv,*tt,*cc,ymin=0;
178
179 vv=(double*)malloc(nnodes*sizeof(double));
180 tt=(double*)malloc(nnodes*sizeof(double));
181 cc=(double*)malloc((nnodes+2)*sizeof(double));
182 info=SDPConeComputeXV(sdpcone,0,&derror);
183 for (i=0;i<nnodes;i++){
184 for (j=0;j<nnodes;j++){dd = (( rand())/scal - .5); vv[j] = tan(3.1415926*dd);}
185 info=SDPConeXVMultiply(sdpcone,0,vv,tt,nnodes);
186 for (j=0;j<nnodes;j++){if (tt[j]<0) tt[j]=-1; else tt[j]=1;}
187 for (j=0;j<nnodes+2;j++){cc[j]=0;}
188 info=SDPConeAddXVAV(sdpcone,0,tt,nnodes,cc,nnodes+2);
189 if (cc[0]<ymin) ymin=cc[0];
190 }
191 printf("Best integer solution: %4.2f\n",ymin);
192 free(vv); free(tt); free(cc);
193
194 return(0);
195 }
196
TCheckArgs(DSDP dsdp,int nargs,char * runargs[])197 static int TCheckArgs(DSDP dsdp,int nargs,char *runargs[]){
198
199 int kk, info;
200
201 info=DSDPSetOptions(dsdp,runargs,nargs);
202 for (kk=1; kk<nargs-1; kk++){
203 if (strncmp(runargs[kk],"-dloginfo",8)==0){
204 info=DSDPLogInfoAllow(DSDP_TRUE,0);
205 } else if (strncmp(runargs[kk],"-params",7)==0){
206 info=DSDPReadOptions(dsdp,runargs[kk+1]);
207 } else if (strncmp(runargs[kk],"-help",7)==0){
208 printf("%s\n",help);
209 }
210 }
211
212 return 0;
213 }
214
215
216 #define BUFFERSIZ 100
217
218 #undef __FUNCT__
219 #define __FUNCT__ "ParseGraphline"
ParseGraphline(char * thisline,int * row,int * col,double * value,int * gotem)220 static int ParseGraphline(char * thisline,int *row,int *col,double *value,
221 int *gotem){
222
223 int temp;
224 int rtmp,coltmp;
225
226 rtmp=-1, coltmp=-1, *value=0.0;
227 temp=sscanf(thisline,"%d %d %lf",&rtmp,&coltmp,value);
228 if (temp==3 && rtmp>0 && coltmp>0) *gotem=3;
229 else if (temp==2 && rtmp>0 && coltmp>0){ *value = 1.0; *gotem=3;}
230 else *gotem=0;
231 *row=rtmp-1; *col=coltmp-1;
232
233 return(0);
234 }
235
236
237 #undef __FUNCT__
238 #define __FUNCT__ "ReadGraph"
ReadGraph(char * filename,int * nnodes,int * nedges,int ** n1,int ** n2,double ** wght)239 int ReadGraph(char* filename,int *nnodes, int *nedges,
240 int**n1, int ** n2, double **wght){
241
242 FILE*fp;
243 char thisline[BUFFERSIZ]="*";
244 int i,k=0,line=0,nodes,edges,gotone=3;
245 int *node1,*node2;
246 double *weight;
247 int info,row,col;
248 double value;
249
250 fp=fopen(filename,"r");
251 if (!fp){printf("Cannot open file %s !",filename);return(1);}
252
253 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){
254 fgets(thisline,BUFFERSIZ,fp); line++;
255 }
256
257 if (sscanf(thisline,"%d %d",&nodes, &edges)!=2){
258 printf("First line must contain the number of nodes and number of edges\n");
259 return 1;
260 }
261
262 node1=(int*)malloc(edges*sizeof(int));
263 node2=(int*)malloc(edges*sizeof(int));
264 weight=(double*)malloc(edges*sizeof(double));
265
266 for (i=0; i<edges; i++){ node1[i]=0;node2[i]=0;weight[i]=0.0;}
267
268 while(!feof(fp) && gotone){
269 thisline[0]='\0';
270 fgets(thisline,BUFFERSIZ,fp); line++;
271 info = ParseGraphline(thisline,&row,&col,&value,&gotone);
272 if (gotone && value!=0.0 && k<edges &&
273 col < nodes && row < nodes && col >= 0 && row >= 0){
274 if (row<col){info=row;row=col;col=info;}
275 if (row == col){}
276 else {
277 node1[k]=row; node2[k]=col;
278 weight[k]=value; k++;
279 }
280 } else if (gotone && k>=edges) {
281 printf("To many edges in file.\nLine %d, %s\n",line,thisline);
282 return 1;
283 } else if (gotone&&(col >= nodes || row >= nodes || col < 0 || row < 0)){
284 printf("Invalid node number.\nLine %d, %s\n",line,thisline);
285 return 1;
286 }
287 }
288 *nnodes=nodes; *nedges=edges;
289 *n1=node1; *n2=node2; *wght=weight;
290 return 0;
291 }
292