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 DSDPSetup(DSDP dsdp)
Set up data structures in the solver and the cones associated with it.
Definition: dsdpsetup.c:193
int DSDPDestroy(DSDP dsdp)
Free the internal data structures of the solver and the cones associated with it.
Definition: dsdpsetup.c:496
int DSDPSetOptions(DSDP dsdp, char *runargs[], int nargs)
Read command line arguments to set options in DSDP.
int DSDPSetDualObjective(DSDP dsdp, int i, double bi)
Set the objective vector b in (D).
Definition: dsdpsetdata.c:25
int DSDPGetY(DSDP dsdp, double y[], int m)
Copies the variables y into an array.
Definition: dsdpsetdata.c:100
int DSDPSetStandardMonitor(DSDP dsdp, int k)
Print at every kth iteration.
Definition: dsdpprintout.c:153
int DSDPCreate(int m, DSDP *dsdpnew)
Create a DSDP solver. FIRST DSDP routine!
Definition: dsdpsetup.c:30
int DSDPSolve(DSDP dsdp)
Apply DSDP to the problem.
Definition: dsdpsetup.c:343
int DSDPComputeX(DSDP dsdp)
Compute the X variables.
Definition: dsdpx.c:55
int DSDPSetGapTolerance(DSDP dsdp, double gaptol)
Terminate the solver when the relative duality gap is less than this tolerance.
Definition: dsdpconverge.c:110
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 DSDPStopReason(DSDP dsdp, DSDPTerminationReason *reason)
Copy the reason why the solver terminated.
Definition: dsdpsetdata.c:582
int DSDPSetR0(DSDP dsdp, double res)
Set an initial value for the variable r in (DD)
Definition: dsdpsetdata.c:311
int DSDPSetY0(DSDP dsdp, int i, double yi0)
Set the initial values of variables y in (D).
Definition: dsdpsetdata.c:77
int DSDPSetZBar(DSDP dsdp, double ppobj)
Set an upper bound on the objective value at the solution.
Definition: dsdpsetdata.c:283
int DSDPSetPotentialParameter(DSDP dsdp, double rho)
Set the potential parameter.
Definition: dsdpsetdata.c:765
int DSDPReuseMatrix(DSDP dsdp, int rm)
Reuse the Hessian of the barrier function multiple times at each DSDP iteration.
Definition: dsdpsetdata.c:905
int DSDPReadOptions(DSDP dsdp, char filename[])
Read DSDP parameters from a file.
int MaxCut(int, int, int[], int[], double[])
Formulate and solve the SDP relaxation of the Maximum Cut problem.
Definition: maxcut.c:51
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
int SDPConeViewDataMatrix(SDPCone sdpcone, int blockj, int vari)
Print a data matrix to the screen.
Definition: dsdpadddata.c:205
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
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 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 SDPConeSetBlockSize(SDPCone sdpcone, int blockj, int n)
Set the dimension of one block in the semidefinite cone.
Definition: dsdpadddata.c:535
int SDPConeComputeXV(SDPCone sdpcone, int blockj, int *derror)
Compute a factor V such that .
Definition: sdpcone.c:325
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 SDPConeAddXVAV(SDPCone sdpcone, int blockj, double vin[], int n, double sum[], int mm)
Compute for i = 0 through m.
Definition: sdpcone.c:292
Internal structures for the DSDP solver.
Definition: dsdp.h:65
Internal structure for semidefinite cone.
Definition: dsdpsdp.h:80