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>
12char help[]="\
13DSDP Usage: maxcut <graph filename> \n\
14 -gaptol <relative duality gap: default is 0.001> \n\
15 -maxit <maximum iterates: default is 200> \n";
16
17static int ReadGraph(char*,int *, int *,int**, int **, double **);
18static int TCheckArgs(DSDP,int,char **);
19static int ParseGraphline(char*,int*,int*,double*,int*);
21int MaxCut(int,int, int[], int[], double[]);
22
23
24int 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
51int 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;
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
175int 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
197static 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"
220static 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"
239int 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}
The API to DSDP for those applications using DSDP as a subroutine library.
DSDPTerminationReason
There are many reasons to terminate the solver.
@ DSDP_INFEASIBLE_START
@ DSDP_TRUE
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
Internal structures for the DSDP solver.
Definition dsdp.h:65
Internal structure for semidefinite cone.
Definition dsdpsdp.h:80