DSDP
|
00001 #include "dsdp5.h" 00002 #include <string.h> 00003 #include <math.h> 00004 #include <stdio.h> 00005 #include <stdlib.h> 00006 00012 char help2[]="\nA positive semidefinite relaxation of the\n\ 00013 graph k coloring problem can be rewritten as\n\n\ 00014 Find X>=0 \n\ 00015 such that X_ij <= 1 - 1/(k-1) for all edges (i,j).\n\ 00016 "; 00017 00018 char help[]="DSDP Usage: color <graph filename> "; 00019 00020 static int ReadGraph(char*,int *, int *,int**, int **, double **); 00021 static int ParseGraphline(char*,int*,int*,double*,int*); 00022 static int RandomizedColor(DSDP, SDPCone, int, int[], int[], int); 00023 int MinColoring(int argc,char *argv[]); 00024 00025 00026 int main(int argc,char *argv[]){ 00027 int info; 00028 info=MinColoring(argc,argv); 00029 return 0; 00030 } 00031 00039 int MinColoring(int argc,char *argv[]){ 00040 00041 int i,kk,vari,info; 00042 int *node1,*node2,nedges,nnodes; 00043 int *iptr1,*iptr2; 00044 double *weight,*yy1,*yy2,bb; 00045 DSDPTerminationReason reason; 00046 SDPCone sdpcone; 00047 BCone bcone; 00048 DSDP dsdp; 00049 00050 if (argc<2){ printf("%s\n%s",help2,help); return(1); } 00051 00052 info = ReadGraph(argv[1],&nnodes,&nedges,&node1,&node2,&weight); 00053 if (info){ printf("Problem reading file\n"); return 1; } 00054 00055 info = DSDPCreate(nnodes+nedges,&dsdp); 00056 00057 info = DSDPCreateSDPCone(dsdp,1,&sdpcone); 00058 info = SDPConeSetBlockSize(sdpcone,0,nnodes); 00059 info = SDPConeSetSparsity(sdpcone,0,nnodes+nedges+1); 00060 00061 info = DSDPCreateBCone(dsdp,&bcone); 00062 00063 if (info){ printf("Out of memory\n"); return 1; } 00064 00065 00066 /* Formulate the problem from the data */ 00067 /* Create data matrices */ 00068 00069 /* Diagonal elements of X(i,i) must equal 1.0 */ 00070 iptr1=(int*)malloc(nnodes*sizeof(int)); 00071 yy1=(double*)malloc(nnodes*sizeof(double)); 00072 for (i=0;i<nnodes;i++){ 00073 iptr1[i]=(i+2)*(i+1)/2-1; 00074 yy1[i]=1.0; 00075 } 00076 for (i=0;i<nnodes;i++){ 00077 info=SDPConeSetASparseVecMat(sdpcone,0,i+1,nnodes,1.0,0,iptr1+i,yy1+i,1); 00078 if (info) printf("ERROR 1: %d \n",i); 00079 info=DSDPSetDualObjective(dsdp,i+1,1.0); 00080 } 00081 00082 /* For each nonzero element (i,j) of the matrix, X(i,j) must be less than 1 - 1/nnodes */ 00083 bb=2-2.0/nnodes; 00084 iptr2=(int*)malloc(nedges*sizeof(int)); 00085 yy2=(double*)malloc(nedges*sizeof(double)); 00086 for (i=0;i<nedges;i++){ 00087 iptr2[i]=(node1[i])*(node1[i]+1)/2+node2[i]; 00088 yy2[i]=1.0; 00089 } 00090 info = BConeAllocateBounds(bcone,nedges); 00091 for (i=0; i<nedges; i++){ 00092 vari=nnodes+i+1; 00093 info = SDPConeSetSparseVecMat(sdpcone,0,vari,nnodes,0,iptr2+i,yy2+i,1); 00094 if (info) printf("ERROR 2: %d %d \n",i,vari); 00095 info = BConeSetPSlackVariable(bcone,vari); 00096 if (info) printf("ERROR 3: %d %d \n",i,vari); 00097 info = DSDPSetDualObjective(dsdp,vari,bb); 00098 } 00099 00100 00101 /* Get read to go */ 00102 info=DSDPSetPotentialParameter(dsdp,5); 00103 00104 for (kk=1; kk<argc-1; kk++){ 00105 if (strncmp(argv[kk],"-dloginfo",8)==0){ 00106 info=DSDPLogInfoAllow(DSDP_TRUE,0); 00107 } else if (strncmp(argv[kk],"-params",7)==0){ 00108 info=DSDPReadOptions(dsdp,argv[kk+1]); 00109 } else if (strncmp(argv[kk],"-help",7)==0){ 00110 printf("%s\n",help); 00111 } 00112 } 00113 info=DSDPSetOptions(dsdp,argv,argc); 00114 00115 if (info){ printf("Out of memory\n"); return 1; } 00116 info = DSDPSetStandardMonitor(dsdp,1); 00117 00118 info = DSDPSetup(dsdp); 00119 if (info){ printf("Out of memory\n"); return 1; } 00120 00121 info = DSDPSolve(dsdp); 00122 if (info){ printf("Numerical error\n"); return 1; } 00123 info = DSDPStopReason(dsdp,&reason); 00124 00125 if (reason!=DSDP_INFEASIBLE_START){ /* Randomized solution strategy */ 00126 info=RandomizedColor(dsdp, sdpcone, nnodes, node1, node2, nedges); 00127 } 00128 00129 info = DSDPDestroy(dsdp); 00130 00131 free(node1);free(node2);free(weight); 00132 free(iptr1); 00133 free(yy1); 00134 free(iptr2); 00135 free(yy2); 00136 00137 return 0; 00138 } 00139 00140 static int GetXRow(double xmat[],double xrow[],int row,int n){ 00141 int i,i1=row*(row+1)/2; 00142 for (i=0;i<row;i++){xrow[i]=xmat[i1+i];} 00143 for (i=row;i<n;i++){xrow[i]=xmat[i1+row];i1+=i+1;} 00144 return 0; 00145 } 00146 00147 typedef struct { 00148 int index;double val; 00149 } orderVec; 00150 00151 static int cut_comp(const void *e1,const void *e2){ /* Used in qsort routine */ 00152 double d1=((orderVec*)e1)->val, d2=((orderVec*)e2)->val; 00153 if (d1<d2) return (1); 00154 else if (d1>d2) return (-1); 00155 return(0); 00156 } 00157 00158 static int RemoveNode(int node, int node1[], int node2[], int *nedges){ 00159 int i,nnedges=*nedges; 00160 for (i=0;i<nnedges;i++){ 00161 if (node1[i]==node || node2[i]==node){ 00162 node1[i]=node1[nnedges-1]; 00163 node2[i]=node2[nnedges-1]; 00164 nnedges--; 00165 if (i < nnedges) i--; 00166 } 00167 } 00168 *nedges=nnedges; 00169 return 0; 00170 } 00171 00172 static int Connected(int n1, int n2, int node1[], int node2[], int nedges){ 00173 int i; 00174 if (n1==n2) return 1; 00175 for (i=0;i<nedges;i++){ 00176 if (node1[i]==n1 && node2[i]==n2){ return 1;} 00177 if (node1[i]==n2 && node2[i]==n1){ return 1;} 00178 } 00179 return 0; 00180 } 00181 00182 static int HighDegree(int node1[], int node2[], int nedges, int iwork[], int nnodes){ 00183 int i,nmax=0,maxdegree=-1; 00184 for (i=0;i<nnodes;i++) iwork[i]=0; 00185 for (i=0;i<nedges;i++){ 00186 iwork[node1[i]]++; iwork[node2[i]]++; 00187 } 00188 for (i=0;i<nnodes;i++){ if (iwork[i]>maxdegree){nmax=i; maxdegree=iwork[i];} } 00189 return nmax; 00190 } 00191 00192 static int First(int coloring[], int nnodes){ 00193 int i,nmax=nnodes; 00194 for (i=0;i<nnodes;i++){ 00195 if (coloring[i]==0){ 00196 nmax=i; return nmax; 00197 } 00198 } 00199 return -1; 00200 } 00201 00202 static int RandomizedColor(DSDP dsdp, SDPCone sdpcone, int nnodes, int node1[], int node2[], int nedges){ 00203 int i,j,nodek,nn,info,flag,coloring=0,maxvertex; 00204 int *degree,*color,*cgroup,ngsize,uncolored=nnodes; 00205 int tnedges=nedges; 00206 double *xrow,*xptr; 00207 orderVec *vorder; 00208 00209 xrow=(double*)malloc(nnodes*sizeof(double)); 00210 color=(int*)malloc(nnodes*sizeof(int)); 00211 cgroup=(int*)malloc(nnodes*sizeof(int)); 00212 degree=(int*)malloc(nnodes*sizeof(int)); 00213 vorder=(orderVec*)malloc(nnodes*sizeof(orderVec)); 00214 00215 for (i=0;i<nnodes;i++){ color[i]=0;} 00216 info=DSDPComputeX(dsdp); 00217 info=SDPConeGetXArray(sdpcone,0,&xptr,&nn); 00218 00219 while (uncolored>0){ 00220 00221 coloring++; 00222 00223 maxvertex=First(color,nnodes); 00224 maxvertex=HighDegree(node1,node2,tnedges,degree,nnodes); 00225 00226 cgroup[0]=maxvertex;ngsize=1; 00227 00228 info=GetXRow(xptr,xrow,maxvertex,nnodes); 00229 00230 for (i=0;i<nnodes;i++){vorder[i].index=i; vorder[i].val = xrow[i];} 00231 qsort( (void*)vorder, nnodes, sizeof(orderVec), cut_comp); 00232 00233 for (i=0;i<nnodes;i++){ 00234 nodek=vorder[i].index; 00235 if (color[nodek]==0){ 00236 for (flag=0,j=0;j<ngsize;j++){ 00237 if (Connected(nodek,cgroup[j],node1,node2,tnedges) ){flag=1;} 00238 } 00239 if (flag==0){ cgroup[ngsize]=nodek; ngsize++; } 00240 } 00241 } 00242 for (i=0;i<ngsize;i++){ 00243 color[cgroup[i]]=coloring; uncolored--; 00244 RemoveNode(cgroup[i],node1,node2,&tnedges); 00245 } 00246 } 00247 printf("\nCOLORS: %d\n",coloring); 00248 free(xrow); 00249 free(color); 00250 free(cgroup); 00251 free(degree); 00252 free(vorder); 00253 return(0); 00254 } 00255 00256 00257 #define BUFFERSIZ 100 00258 00259 #undef __FUNCT__ 00260 #define __FUNCT__ "ParseGraphline" 00261 static int ParseGraphline(char * thisline,int *row,int *col,double *value, 00262 int *gotem){ 00263 00264 int temp; 00265 int rtmp,coltmp; 00266 00267 rtmp=-1, coltmp=-1, *value=0.0; 00268 temp=sscanf(thisline,"%d %d %lf",&rtmp,&coltmp,value); 00269 if (temp==3 && rtmp>0 && coltmp>0) *gotem=3; 00270 else if (temp==2 && rtmp>0 && coltmp>0){ *value = 1.0; *gotem=3;} 00271 else *gotem=0; 00272 *row=rtmp-1; *col=coltmp-1; 00273 00274 return(0); 00275 } 00276 00277 00278 #undef __FUNCT__ 00279 #define __FUNCT__ "ReadGraph" 00280 int ReadGraph(char* filename,int *nnodes, int *nedges, 00281 int**n1, int ** n2, double **wght){ 00282 00283 FILE*fp; 00284 char thisline[BUFFERSIZ]="*"; 00285 int i,k=0,line=0,nodes,edges,gotone=3; 00286 int *node1,*node2; 00287 double *weight; 00288 int info,row,col; 00289 double value; 00290 00291 fp=fopen(filename,"r"); 00292 if (!fp){printf("Cannot open file %s !",filename);return(1);} 00293 00294 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){ 00295 fgets(thisline,BUFFERSIZ,fp); line++; 00296 } 00297 00298 if (sscanf(thisline,"%d %d",&nodes, &edges)!=2){ 00299 printf("First line must contain the number of nodes and number of edges\n"); 00300 return 1; 00301 } 00302 00303 node1=(int*)malloc(edges*sizeof(int)); 00304 node2=(int*)malloc(edges*sizeof(int)); 00305 weight=(double*)malloc(edges*sizeof(double)); 00306 00307 for (i=0; i<edges; i++){ node1[i]=0;node2[i]=0;weight[i]=0.0;} 00308 00309 while(!feof(fp) && gotone){ 00310 thisline[0]='\0'; 00311 fgets(thisline,BUFFERSIZ,fp); line++; 00312 info = ParseGraphline(thisline,&row,&col,&value,&gotone); 00313 if (gotone && value!=0.0 && k<edges && 00314 col < nodes && row < nodes && col >= 0 && row >= 0){ 00315 if (row<col){info=row;row=col;col=info;} 00316 if (row == col){} 00317 else { 00318 node1[k]=row; node2[k]=col; 00319 weight[k]=value; k++; 00320 } 00321 } else if (gotone && k>=edges) { 00322 printf("To many edges in file.\nLine %d, %s\n",line,thisline); 00323 return 1; 00324 } else if (gotone&&(col >= nodes || row >= nodes || col < 0 || row < 0)){ 00325 printf("Invalid node number.\nLine %d, %s\n",line,thisline); 00326 return 1; 00327 } 00328 } 00329 *nnodes=nodes; *nedges=edges; 00330 *n1=node1; *n2=node2; *wght=weight; 00331 return 0; 00332 }