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>
11
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.
struct SDPCone_C * SDPCone
The SDPCone object points to blocks of data that specify semidefinite matrix inequalities.
Definition dsdp5.h:26
@ DSDP_TRUE
struct DSDP_C * DSDP
An implementation of the dual-scaling algorithm for semidefinite programming.
DSDPTerminationReason
There are many reasons to terminate the solver.
@ DSDP_INFEASIBLE_START
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.
int DSDPSetStandardMonitor(DSDP dsdp, int k)
Print at every kth iteration.
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.
int DSDPSetPNormTolerance(DSDP dsdp, double ptol)
Terminate the solver when the relative duality gap is suffiently small and the PNorm is less than thi...
int DSDPStopReason(DSDP dsdp, DSDPTerminationReason *reason)
Copy the reason why the solver terminated.
int DSDPSetR0(DSDP dsdp, double res)
Set an initial value for the variable r in (DD)
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.
int DSDPSetPotentialParameter(DSDP dsdp, double rho)
Set the potential parameter.
int DSDPReuseMatrix(DSDP dsdp, int rm)
Reuse the Hessian of the barrier function multiple times at each DSDP iteration.
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.
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.
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.
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