DSDP
color.c
Go to the documentation of this file.
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 }