DSDP
|
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 /* Formulate the problem from the data */ 00069 /* 00070 Diagonal elements equal 1.0 00071 Create Constraint matrix A_i for i=1, ..., nnodes. 00072 that has a single nonzero element. 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 /* C matrix is the Laplacian of the adjacency matrix */ 00091 /* Also compute a feasible initial point y such that S>=0 */ 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 { /* Equivalent */ 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 /* Initial Point */ 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 /* Get read to go */ 00126 info=DSDPSetGapTolerance(dsdp,0.001); 00127 info=DSDPSetPotentialParameter(dsdp,5); 00128 info=DSDPReuseMatrix(dsdp,0); 00129 info=DSDPSetPNormTolerance(dsdp,1.0); 00130 /* 00131 info = TCheckArgs(dsdp,argc,argv); 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){ /* Randomized solution strategy */ 00145 info=MaxCutRandomized(sdpcone,nnodes); 00146 if (0==1){ /* Look at the solution */ 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 } /* main */ 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 }