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 }
int DSDPCreate(int m, DSDP *dsdpnew)
Create a DSDP solver. FIRST DSDP routine!
Definition: dsdpsetup.c:30
int DSDPDestroy(DSDP dsdp)
Free the internal data structures of the solver and the cones associated with it. ...
Definition: dsdpsetup.c:496
int DSDPReadOptions(DSDP dsdp, char filename[])
Read DSDP parameters from a file.
int DSDPSetOptions(DSDP dsdp, char *runargs[], int nargs)
Read command line arguments to set options in DSDP.
int DSDPSetZBar(DSDP dsdp, double ppobj)
Set an upper bound on the objective value at the solution.
Definition: dsdpsetdata.c:283
int DSDPSetDualObjective(DSDP dsdp, int i, double bi)
Set the objective vector b in (D).
Definition: dsdpsetdata.c:25
Internal structures for the DSDP solver.
Definition: dsdp.h:65
DSDPTerminationReason
There are many reasons to terminate the solver.
int DSDPSetGapTolerance(DSDP dsdp, double gaptol)
Terminate the solver when the relative duality gap is less than this tolerance.
Definition: dsdpconverge.c:110
int SDPConeGetXArray(SDPCone sdpcone, int blockj, double *xx[], int *nn)
After applying the solver, set a pointer to the array in the object with the solution X...
Definition: dsdpadddata.c:328
The API to DSDP for those applications using DSDP as a subroutine library.
int SDPConeAddXVAV(SDPCone sdpcone, int blockj, double vin[], int n, double sum[], int mm)
Compute for i = 0 through m.
Definition: sdpcone.c:292
int SDPConeViewDataMatrix(SDPCone sdpcone, int blockj, int vari)
Print a data matrix to the screen.
Definition: dsdpadddata.c:205
int DSDPSetup(DSDP dsdp)
Set up data structures in the solver and the cones associated with it.
Definition: dsdpsetup.c:193
int DSDPReuseMatrix(DSDP dsdp, int rm)
Reuse the Hessian of the barrier function multiple times at each DSDP iteration.
Definition: dsdpsetdata.c:905
int MaxCut(int, int, int[], int[], double[])
Formulate and solve the SDP relaxation of the Maximum Cut problem.
Definition: maxcut.c:51
int SDPConeAddASparseVecMat(SDPCone sdpcone, int blockj, int vari, int n, double alpha, int ishift, const int ind[], const double val[], int nnz)
Add data matrix in a sparse format.
int DSDPSetPotentialParameter(DSDP dsdp, double rho)
Set the potential parameter.
Definition: dsdpsetdata.c:765
int DSDPStopReason(DSDP dsdp, DSDPTerminationReason *reason)
Copy the reason why the solver terminated.
Definition: dsdpsetdata.c:582
int SDPConeComputeXV(SDPCone sdpcone, int blockj, int *derror)
Compute a factor V such that .
Definition: sdpcone.c:325
int DSDPSetY0(DSDP dsdp, int i, double yi0)
Set the initial values of variables y in (D).
Definition: dsdpsetdata.c:77
int DSDPSetPNormTolerance(DSDP dsdp, double ptol)
Terminate the solver when the relative duality gap is suffiently small and the PNorm is less than thi...
Definition: dsdpconverge.c:158
int SDPConeXVMultiply(SDPCone sdpcone, int blockj, double vin[], double vout[], int n)
Multiply an array by a factor V such that .
Definition: sdpcone.c:251
int DSDPSolve(DSDP dsdp)
Apply DSDP to the problem.
Definition: dsdpsetup.c:343
int SDPConeSetBlockSize(SDPCone sdpcone, int blockj, int n)
Set the dimension of one block in the semidefinite cone.
Definition: dsdpadddata.c:535
Internal structure for semidefinite cone.
Definition: dsdpsdp.h:80
int DSDPGetY(DSDP dsdp, double y[], int m)
Copies the variables y into an array.
Definition: dsdpsetdata.c:100
int DSDPComputeX(DSDP dsdp)
Compute the X variables.
Definition: dsdpx.c:55
int SDPConeSetASparseVecMat(SDPCone sdpcone, int blockj, int vari, int n, double alpha, int ishift, const int ind[], const double val[], int nnz)
Set data matrix in a sparse format.
int DSDPSetR0(DSDP dsdp, double res)
Set an initial value for the variable r in (DD)
Definition: dsdpsetdata.c:311
int DSDPSetStandardMonitor(DSDP dsdp, int k)
Print at every kth iteration.
Definition: dsdpprintout.c:153
int MaxCutRandomized(SDPCone, int)
Apply the Goemens and Williamson randomized cut algorithm to the SDP relaxation of the max-cut proble...
Definition: maxcut.c:175