00001 #include "dsdp5.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <stdio.h>
00005 #include <stdlib.h>
00012 char help[]="\
00013 DSDP Usage: maxcut <graph filename> \n\
00014 -gaptol <relative duality gap: default is 0.001> \n\
00015 -maxit <maximum iterates: default is 200> \n";
00016
00017 static int ReadGraph(char*,int *, int *,int**, int **, double **);
00018 static int TCheckArgs(DSDP,int,char **);
00019 static int ParseGraphline(char*,int*,int*,double*,int*);
00020 int MaxCutRandomized(SDPCone,int);
00021 int MaxCut(int,int, int[], int[], double[]);
00022
00023
00024 int main(int argc,char *argv[]){
00025 int info;
00026 int *node1,*node2,nedges,nnodes;
00027 double *weight;
00028
00029 if (argc<2){ printf("%s",help); return(1); }
00030
00031 info = ReadGraph(argv[1],&nnodes,&nedges,&node1,&node2,&weight);
00032 if (info){ printf("Problem reading file\n"); return 1; }
00033
00034 MaxCut(nnodes,nedges,node1,node2,weight);
00035
00036 free(node1);free(node2);free(weight);
00037 return 0;
00038 }
00039
00051 int MaxCut(int nnodes,int nedges, int node1[], int node2[], double weight[]){
00052
00053 int i,info;
00054 int *indd,*iptr;
00055 double *yy,*val,*diag,tval=0;
00056 DSDPTerminationReason reason;
00057 SDPCone sdpcone;
00058 DSDP dsdp;
00059
00060 info = DSDPCreate(nnodes,&dsdp);
00061 info = DSDPCreateSDPCone(dsdp,1,&sdpcone);
00062
00063 if (info){ printf("Out of memory\n"); return 1; }
00064
00065 info = SDPConeSetBlockSize(sdpcone,0,nnodes);
00066
00067
00068
00069
00070
00071
00072
00073
00074 diag=(double*)malloc(nnodes*sizeof(double));
00075 iptr=(int*)malloc(nnodes*sizeof(int));
00076 for (i=0;i<nnodes;i++){
00077 iptr[i]=i*(i+1)/2+i;
00078 diag[i]=1.0;
00079 }
00080
00081 for (i=0;i<nnodes;i++){
00082 info = DSDPSetDualObjective(dsdp,i+1,1.0);
00083 info = SDPConeSetASparseVecMat(sdpcone,0,i+1,nnodes,1.0,0,iptr+i,diag+i,1);
00084 if (0==1){
00085 printf("Matrix: %d\n",i+1);
00086 info = SDPConeViewDataMatrix(sdpcone,0,i+1);
00087 }
00088 }
00089
00090
00091
00092 yy=(double*)malloc(nnodes*sizeof(double));
00093 for (i=0;i<nnodes;i++){yy[i]=0.0;}
00094 indd=(int*)malloc((nnodes+nedges)*sizeof(int));
00095 val=(double*)malloc((nnodes+nedges)*sizeof(double));
00096 for (i=0;i<nnodes+nedges;i++){indd[i]=0;}
00097 for (i=0;i<nnodes;i++){indd[nedges+i]=i*(i+1)/2+i;}
00098 for (i=0;i<nnodes+nedges;i++){val[i]=0;}
00099 for (i=0;i<nedges;i++){
00100 indd[i]=(node1[i])*(node1[i]+1)/2 + node2[i];
00101 tval+=fabs(weight[i]);
00102 val[i]=weight[i]/4;
00103 val[nedges+node1[i]]-=weight[i]/4;
00104 val[nedges+node2[i]]-=weight[i]/4;
00105 yy[node1[i]]-= fabs(weight[i]/2);
00106 yy[node2[i]]-= fabs(weight[i]/2);
00107 }
00108
00109 if (0){
00110 info = SDPConeSetASparseVecMat(sdpcone,0,0,nnodes,1.0,0,indd,val,nedges+nnodes);
00111 } else {
00112 info = SDPConeSetASparseVecMat(sdpcone,0,0,nnodes,1.0,0,indd,val,nedges);
00113 info = SDPConeAddASparseVecMat(sdpcone,0,0,nnodes,1.0,0,indd+nedges,val+nedges,nnodes);
00114 }
00115 if (0==1){ info = SDPConeViewDataMatrix(sdpcone,0,0);}
00116
00117
00118 info = DSDPSetR0(dsdp,0.0);
00119 info = DSDPSetZBar(dsdp,10*tval+1.0);
00120 for (i=0;i<nnodes; i++){
00121 info = DSDPSetY0(dsdp,i+1,10*yy[i]);
00122 }
00123 if (info) return info;
00124
00125
00126 info=DSDPSetGapTolerance(dsdp,0.001);
00127 info=DSDPSetPotentialParameter(dsdp,5);
00128 info=DSDPReuseMatrix(dsdp,0);
00129 info=DSDPSetPNormTolerance(dsdp,1.0);
00130
00131
00132
00133
00134 if (info){ printf("Out of memory\n"); return 1; }
00135 info = DSDPSetStandardMonitor(dsdp,1);
00136
00137 info = DSDPSetup(dsdp);
00138 if (info){ printf("Out of memory\n"); return 1; }
00139
00140 info = DSDPSolve(dsdp);
00141 if (info){ printf("Numerical error\n"); return 1; }
00142 info = DSDPStopReason(dsdp,&reason);
00143
00144 if (reason!=DSDP_INFEASIBLE_START){
00145 info=MaxCutRandomized(sdpcone,nnodes);
00146 if (0==1){
00147 int n; double *xx,*y=diag;
00148 info=DSDPGetY(dsdp,y,nnodes);
00149 info=DSDPComputeX(dsdp);DSDPCHKERR(info);
00150 info=SDPConeGetXArray(sdpcone,0,&xx,&n);
00151 }
00152 }
00153 info = DSDPDestroy(dsdp);
00154
00155 free(iptr);
00156 free(yy);
00157 free(indd);
00158 free(val);
00159 free(diag);
00160
00161 return 0;
00162 }
00163
00164
00165
00175 int MaxCutRandomized(SDPCone sdpcone,int nnodes){
00176 int i,j,derror,info;
00177 double dd,scal=RAND_MAX,*vv,*tt,*cc,ymin=0;
00178
00179 vv=(double*)malloc(nnodes*sizeof(double));
00180 tt=(double*)malloc(nnodes*sizeof(double));
00181 cc=(double*)malloc((nnodes+2)*sizeof(double));
00182 info=SDPConeComputeXV(sdpcone,0,&derror);
00183 for (i=0;i<nnodes;i++){
00184 for (j=0;j<nnodes;j++){dd = (( rand())/scal - .5); vv[j] = tan(3.1415926*dd);}
00185 info=SDPConeXVMultiply(sdpcone,0,vv,tt,nnodes);
00186 for (j=0;j<nnodes;j++){if (tt[j]<0) tt[j]=-1; else tt[j]=1;}
00187 for (j=0;j<nnodes+2;j++){cc[j]=0;}
00188 info=SDPConeAddXVAV(sdpcone,0,tt,nnodes,cc,nnodes+2);
00189 if (cc[0]<ymin) ymin=cc[0];
00190 }
00191 printf("Best integer solution: %4.2f\n",ymin);
00192 free(vv); free(tt); free(cc);
00193
00194 return(0);
00195 }
00196
00197 static int TCheckArgs(DSDP dsdp,int nargs,char *runargs[]){
00198
00199 int kk, info;
00200
00201 info=DSDPSetOptions(dsdp,runargs,nargs);
00202 for (kk=1; kk<nargs-1; kk++){
00203 if (strncmp(runargs[kk],"-dloginfo",8)==0){
00204 info=DSDPLogInfoAllow(DSDP_TRUE,0);
00205 } else if (strncmp(runargs[kk],"-params",7)==0){
00206 info=DSDPReadOptions(dsdp,runargs[kk+1]);
00207 } else if (strncmp(runargs[kk],"-help",7)==0){
00208 printf("%s\n",help);
00209 }
00210 }
00211
00212 return 0;
00213 }
00214
00215
00216 #define BUFFERSIZ 100
00217
00218 #undef __FUNCT__
00219 #define __FUNCT__ "ParseGraphline"
00220 static int ParseGraphline(char * thisline,int *row,int *col,double *value,
00221 int *gotem){
00222
00223 int temp;
00224 int rtmp,coltmp;
00225
00226 rtmp=-1, coltmp=-1, *value=0.0;
00227 temp=sscanf(thisline,"%d %d %lf",&rtmp,&coltmp,value);
00228 if (temp==3 && rtmp>0 && coltmp>0) *gotem=3;
00229 else if (temp==2 && rtmp>0 && coltmp>0){ *value = 1.0; *gotem=3;}
00230 else *gotem=0;
00231 *row=rtmp-1; *col=coltmp-1;
00232
00233 return(0);
00234 }
00235
00236
00237 #undef __FUNCT__
00238 #define __FUNCT__ "ReadGraph"
00239 int ReadGraph(char* filename,int *nnodes, int *nedges,
00240 int**n1, int ** n2, double **wght){
00241
00242 FILE*fp;
00243 char thisline[BUFFERSIZ]="*";
00244 int i,k=0,line=0,nodes,edges,gotone=3;
00245 int *node1,*node2;
00246 double *weight;
00247 int info,row,col;
00248 double value;
00249
00250 fp=fopen(filename,"r");
00251 if (!fp){printf("Cannot open file %s !",filename);return(1);}
00252
00253 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){
00254 fgets(thisline,BUFFERSIZ,fp); line++;
00255 }
00256
00257 if (sscanf(thisline,"%d %d",&nodes, &edges)!=2){
00258 printf("First line must contain the number of nodes and number of edges\n");
00259 return 1;
00260 }
00261
00262 node1=(int*)malloc(edges*sizeof(int));
00263 node2=(int*)malloc(edges*sizeof(int));
00264 weight=(double*)malloc(edges*sizeof(double));
00265
00266 for (i=0; i<edges; i++){ node1[i]=0;node2[i]=0;weight[i]=0.0;}
00267
00268 while(!feof(fp) && gotone){
00269 thisline[0]='\0';
00270 fgets(thisline,BUFFERSIZ,fp); line++;
00271 info = ParseGraphline(thisline,&row,&col,&value,&gotone);
00272 if (gotone && value!=0.0 && k<edges &&
00273 col < nodes && row < nodes && col >= 0 && row >= 0){
00274 if (row<col){info=row;row=col;col=info;}
00275 if (row == col){}
00276 else {
00277 node1[k]=row; node2[k]=col;
00278 weight[k]=value; k++;
00279 }
00280 } else if (gotone && k>=edges) {
00281 printf("To many edges in file.\nLine %d, %s\n",line,thisline);
00282 return 1;
00283 } else if (gotone&&(col >= nodes || row >= nodes || col < 0 || row < 0)){
00284 printf("Invalid node number.\nLine %d, %s\n",line,thisline);
00285 return 1;
00286 }
00287 }
00288 *nnodes=nodes; *nedges=edges;
00289 *n1=node1; *n2=node2; *wght=weight;
00290 return 0;
00291 }