DSDP
maxcut.c
Go to the documentation of this file.
1 #include "dsdp5.h"
2 #include <string.h>
3 #include <math.h>
4 #include <stdio.h>
5 #include <stdlib.h>
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 
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 
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 
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 
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"
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"
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 }