theta.c

Go to the documentation of this file.
00001 #include "dsdp5.h"
00002 #include <stdio.h>
00003 #include <string.h>
00004 #include <math.h>
00005 #include <stdlib.h>
00006 
00011 char help[]="\n\
00012 Compute the Lovasz theta number for a graph.  This number is an upper bound for \n\
00013 the maximum clique of a graph, a lower bound for the mimimal graph coloring, and serves \n\
00014 as a bound for several other combitorial graph problems.  The number is the solution to \n\
00015 a semidfinite program. \n\n\
00016 The file should be the complement of the graph \n\
00017 This file also demonstrates the use of customized data matrices in DSDP.\n\n\
00018 DSDP Usage: theta <graph complement filename> \n";
00019 
00020 typedef struct {
00021   int v1,v2;
00022 } EdgeMat;
00023 
00024 #include "../src/sdp/dsdpdatamat_impl.h"
00025 extern int SDPConeAddDataMatrix(SDPCone,int, int, int, char, struct DSDPDataMat_Ops*, void*);
00026 
00027 
00028 /* These variables and routines are for customized SDP Data Matrices */
00029 static struct  DSDPDataMat_Ops edgematop;
00030 static int EdgeMatOperationsInitialize(struct  DSDPDataMat_Ops*);
00031 
00032 static struct  DSDPDataMat_Ops onematops;
00033 static int OneMatOpsInitialize(struct  DSDPDataMat_Ops*);
00034 
00035 static struct  DSDPDataMat_Ops identitymatops;
00036 static int IdentitymatOperationsInitialize(struct  DSDPDataMat_Ops*);
00037 
00038 static int ReadGraphFromFile(char*,int*, int*, EdgeMat*[]); 
00039 int SetThetaData(DSDP, SDPCone, int, int, EdgeMat[]);
00040 int LovaszTheta(int argc,char *argv[]);
00041 
00042 int main(int argc,char *argv[]){
00043   int info;
00044   info=LovaszTheta(argc,argv);
00045   return 0;
00046 }
00047 
00056 int LovaszTheta(int argc,char *argv[]){
00057 
00058   int info,kk,nedges,nnodes;
00059   EdgeMat *Edges;
00060   DSDP dsdp;
00061   SDPCone  sdpcone;
00062 
00063   if (argc<2){ printf("%s",help); return(1); }
00064 
00065   info = ReadGraphFromFile(argv[1],&nnodes,&nedges,&Edges);
00066   if (info){ printf("Problem reading file\n"); return 1; }
00067 
00068   info = DSDPCreate(nedges+1,&dsdp); 
00069   info = DSDPCreateSDPCone(dsdp,1,&sdpcone);
00070   info = SDPConeSetBlockSize(sdpcone,0,nnodes);
00071   info = SDPConeUsePackedFormat(sdpcone,0);
00072   info = SDPConeSetSparsity(sdpcone,0,nedges+1);
00073   if (info){ printf("Out of memory\n"); return 1; }
00074   info = SetThetaData(dsdp, sdpcone, nnodes, nedges, Edges);
00075   if (info){ printf("Out of memory\n"); return 1; }
00076 
00077   info = DSDPSetGapTolerance(dsdp,0.001);
00078   info = DSDPSetZBar(dsdp,nnodes+1);
00079   info = DSDPReuseMatrix(dsdp,10);
00080 
00081   for (kk=1; kk<argc-1; kk++){
00082     if (strncmp(argv[kk],"-dloginfo",8)==0){
00083       info=DSDPLogInfoAllow(DSDP_TRUE,0);
00084     } else if (strncmp(argv[kk],"-params",7)==0){
00085       info=DSDPReadOptions(dsdp,argv[kk+1]);
00086     } else if (strncmp(argv[kk],"-help",5)==0){
00087       printf("%s\n",help);
00088     } 
00089   }
00090   info=DSDPSetOptions(dsdp,argv,argc);
00091 
00092   if (info){ printf("Out of memory\n"); return 1; }
00093   info = DSDPSetStandardMonitor(dsdp,1);
00094   if (info){ printf("Monitor Problem \n"); return 1; }
00095 
00096   info = DSDPSetup(dsdp);
00097   if (info){ printf("Out of memory\n"); return 1; }
00098   if (0==1){info=SDPConeCheckData(sdpcone);}
00099   
00100   info = DSDPSolve(dsdp); 
00101   if (info){ printf("Numberical error\n"); return 1; }
00102   info=DSDPComputeX(dsdp);DSDPCHKERR(info);
00103   
00104   if (0==1){ /* Look at the solution */
00105     int n; double *xx;
00106     info=SDPConeGetXArray(sdpcone,0,&xx,&n);
00107   }
00108 
00109   info = DSDPDestroy(dsdp);
00110   free(Edges);
00111   
00112   return 0;
00113 } /* main */
00114 
00126 int SetThetaData(DSDP dsdp, SDPCone sdpcone, int nodes, int edges, EdgeMat Edge[]){
00127 
00128   int i,info;
00129 
00130   /* Create data matrices -  these are all custom types */
00131 
00132   /* The c matrix has all elements equal to 1.0 */
00133   info=OneMatOpsInitialize(&onematops);
00134   info=SDPConeAddDataMatrix(sdpcone,0,0,nodes,'P',&onematops,0);
00135 
00136   /* For each edge connecting nodes i and j, X(i,j)=X(j,i)=0 */
00137   info=EdgeMatOperationsInitialize(&edgematop);
00138   for (i=0; i<edges; i++){
00139     info = SDPConeAddDataMatrix(sdpcone,0,i+1,nodes,'P',&edgematop,(void*)&Edge[i]);
00140     info = DSDPSetDualObjective(dsdp,i+1,0.0);
00141     info = DSDPSetY0(dsdp,i+1,0.0);
00142     if (info) return info; 
00143   }
00144 
00145   /* The trace of X must equal 1.0 */
00146   info = IdentitymatOperationsInitialize(&identitymatops);
00147   info = SDPConeAddDataMatrix(sdpcone,0,edges+1,nodes,'P',&identitymatops,0);
00148   info = DSDPSetDualObjective(dsdp,edges+1,1.0);
00149   info = DSDPSetY0(dsdp,edges+1,-10*nodes-1);
00150 
00151   /* The initial point y is feasible and near the central path */
00152   info = DSDPSetR0(dsdp,0);
00153   
00154   return(0);
00155 }
00156 
00157 #define BUFFERSIZ 100
00158 
00159 #undef __FUNCT__  
00160 #define __FUNCT__ "ParseGraphline"
00161 static int ParseGraphline(char * thisline,int *row,int *col,double *value, 
00162                           int *gotem){
00163 
00164   int temp;
00165   int rtmp,coltmp;
00166 
00167   rtmp=-1, coltmp=-1, *value=0.0;
00168   temp=sscanf(thisline,"%d %d %lf",&rtmp,&coltmp,value);
00169   if (temp==3 && rtmp>0 && coltmp>0) *gotem=3;
00170   else if (temp==2 && rtmp>0 && coltmp>0){ *value = 1.0; *gotem=3;}
00171   else *gotem=0;
00172   *row=rtmp-1; *col=coltmp-1;
00173   if (*gotem && (*col < 0 || *row < 0)){
00174     printf("Node Number must be positive.\n, %s\n",thisline);
00175   }
00176   return(0);
00177 }
00178 
00179 #undef __FUNCT__  
00180 #define __FUNCT__ "ReadGraphFromFile"
00181 static int ReadGraphFromFile(char* filename,int *nnodes, int *nedges, EdgeMat**EE){
00182 
00183   FILE*fp;
00184   char thisline[BUFFERSIZ]="*";
00185   int i,k=0,line=0,nodes,edges,gotone=3;
00186   int info,row,col;
00187   double value;
00188   EdgeMat *E;
00189 
00190   fp=fopen(filename,"r");
00191   if (!fp){ printf("Cannot open file %s !",filename); return(1); }
00192 
00193   while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){
00194     fgets(thisline,BUFFERSIZ,fp); line++; }
00195 
00196   if (sscanf(thisline,"%d %d",&nodes, &edges)!=2){
00197     printf("First line must contain the number of nodes and number of edges\n");
00198     return 1;
00199   }
00200 
00201   E=(EdgeMat*)malloc(edges*sizeof(EdgeMat)); *EE=E;
00202   for (i=0; i<edges; i++){ E[i].v1=0; E[i].v2=0; }
00203 
00204   while(!feof(fp) && gotone){
00205     thisline[0]='\0';
00206     fgets(thisline,BUFFERSIZ,fp); line++;
00207     info = ParseGraphline(thisline,&row,&col,&value,&gotone);
00208     if (gotone && k<edges && 
00209         col < nodes && row < nodes && col >= 0 && row >= 0){
00210       if (row > col){i=row;row=col;col=i;}
00211       if (row == col){}
00212       else { E[k].v1=row; E[k].v2=col; k++;}
00213     } else if (gotone &&  k>=edges) {
00214       printf("To many edges in file.\nLine %d, %s\n",line,thisline);
00215       return 1;
00216     } else if (gotone&&(col >= nodes || row >= nodes || col < 0 || row < 0)){
00217       printf("Invalid node number.\nLine %d, %s\n",line,thisline);
00218       return 1;
00219     }
00220   }
00221 
00222   *nnodes=nodes; *nedges=k;
00223                                                  
00224   return 0;
00225 }
00226 
00227 
00228 /* SPecial Matrix where each edge represents a constraint matrix */
00229 static int EdgeMatDestroy(void*);
00230 static int EdgeMatView(void*);
00231 static int EdgeMatVecVec(void*, double[], int, double *);
00232 static int EdgeMatDot(void*, double[], int, int, double *);
00233 static int EdgeMatGetRank(void*, int*, int);
00234 static int EdgeMatFactor(void*);
00235 static int EdgeMatGetEig(void*, int, double*, double[], int, int[], int*);
00236 static int EdgeMatAddRowMultiple(void*, int, double, double[], int);
00237 static int EdgeMatAddMultiple(void*, double, double[], int, int);
00238 static int EdgeMatGetRowNnz(void*, int, int[], int*, int);
00239 
00240 static int EdgeMatDestroy(void* AA){
00241   return 0;
00242 }
00243 static int EdgeMatVecVec(void* A, double x[], int n, double *v){
00244   EdgeMat*E=(EdgeMat*)A;
00245   *v=2*x[E->v1]*x[E->v2];
00246   return 0;
00247 }
00248 static int EdgeMatDot(void* A, double x[], int nn, int n, double *v){
00249   EdgeMat*E=(EdgeMat*)A;
00250   int k=E->v2*(E->v2+1)/2 + E->v1;
00251   *v=2*x[k];
00252   return 0;
00253 }
00254 static int EdgeMatView(void* A){
00255   EdgeMat*E=(EdgeMat*)A;
00256   printf(" Row: %d, Column: %d\n",E->v1,E->v2);
00257   printf(" Row: %d, Column: %d\n",E->v2,E->v1);
00258   return 0;
00259 }
00260 static int EdgeMatFactor(void* A){
00261   return 0;
00262 }
00263 static int EdgeMatGetRank(void *A, int*rank, int n){
00264   *rank=2;
00265   return 0;
00266 }
00267 static int EdgeMatGetEig(void*A, int neig, double *eig, double v[], int n,int spind[], int *nind){
00268   EdgeMat*E=(EdgeMat*)A;
00269   double tt=1.0/sqrt(2.0);
00270   memset((void*)v,0,(n)*sizeof(double)); 
00271   memset((void*)spind,0,(n)*sizeof(int)); 
00272   if (neig==0){
00273     v[E->v1]=tt;v[E->v2]=tt;*eig=1;
00274     spind[0]=E->v1;spind[1]=E->v2; *nind=2;
00275   } else if (neig==1){
00276     v[E->v1]=-tt;v[E->v2]=tt;*eig=-1;
00277     spind[0]=E->v1;spind[1]=E->v2; *nind=2;
00278   } else { *eig=0;*nind=0;}
00279   return 0;
00280 }
00281 static int EdgeMatGetRowNnz(void*A, int nrow, int nz[], int *nnzz, int n){
00282   EdgeMat*E=(EdgeMat*)A;
00283   if (nrow==E->v1){ nz[E->v2]++; *nnzz=1;} 
00284   else if (nrow==E->v2){nz[E->v1]++; *nnzz=1;} 
00285   else {*nnzz=0;}
00286   return 0;
00287 }
00288 static int EdgeMatAddRowMultiple(void*A, int nrow, double dd, double rrow[], int n){
00289   EdgeMat*E=(EdgeMat*)A;
00290   if (nrow==E->v1){ rrow[E->v2]+=dd;} 
00291   else if (nrow==E->v2){rrow[E->v1]+=dd;} 
00292   return 0;
00293 }
00294 static int EdgeMatAddMultiple(void*A, double dd, double vv[], int nn, int n){
00295   EdgeMat*E=(EdgeMat*)A;
00296   int k=E->v2*(E->v2+1)/2 + E->v1;
00297   vv[k]+=dd;
00298   return 0;
00299 }
00300 static int EdgeMatFNorm(void*A, int n, double *fnorm){
00301   *fnorm=2.0;
00302   return 0;
00303 }
00304 static int EdgeMatCountNonzeros(void*A, int *nnz, int n){
00305   *nnz=1;
00306   return 0;
00307 }
00308 static const char *datamatname="THETA EDGE MATRIX";
00309 static int EdgeMatOperationsInitialize(struct  DSDPDataMat_Ops* edgematoperator){
00310   int info;
00311   if (edgematoperator==NULL) return 0;
00312   info=DSDPDataMatOpsInitialize(edgematoperator); if (info){ return info;}
00313   edgematoperator->matfactor1=EdgeMatFactor;
00314   edgematoperator->matgetrank=EdgeMatGetRank;
00315   edgematoperator->matgeteig=EdgeMatGetEig;
00316   edgematoperator->matvecvec=EdgeMatVecVec;
00317   edgematoperator->matrownz=EdgeMatGetRowNnz;
00318   edgematoperator->matdot=EdgeMatDot;
00319   edgematoperator->matfnorm2=EdgeMatFNorm;
00320   edgematoperator->matnnz=EdgeMatCountNonzeros;
00321   edgematoperator->mataddrowmultiple=EdgeMatAddRowMultiple;
00322   edgematoperator->mataddallmultiple=EdgeMatAddMultiple;
00323   edgematoperator->matdestroy=EdgeMatDestroy;
00324   edgematoperator->matview=EdgeMatView;
00325   edgematoperator->matname=datamatname;
00326   edgematoperator->id=25;
00327   return 0;
00328 }
00329 
00330 /* SPecial Matrix where all elements equal negative one */
00331 static int OneMatDestroy(void*);
00332 static int OneMatView(void*);
00333 static int OneMatVecVec(void*, double[], int, double *);
00334 static int OneMatDot(void*, double[], int, int, double *);
00335 static int OneMatGetRank(void*, int*, int);
00336 static int OneMatFactor(void*);
00337 static int OneMatGetEig(void*, int, double*, double[], int, int[], int*);
00338 static int OneMatRowNnz(void*, int, int[], int*, int);
00339 static int OneMatAddRowMultiple(void*, int, double, double[], int);
00340 static int OneMatAddMultiple(void*, double, double[], int,int);
00341 
00342 
00343 static int OneMatFactor(void*A){return 0;}
00344 static int OneMatGetRank(void *A, int *rank, int n){*rank=1;return 0;}
00345 static int OneMatFNorm2(void*AA, int n, double *fnorm2){*fnorm2=1.0*n*n;return 0;}
00346 static int OneMatCountNonzeros(void*AA, int *nnz, int n){*nnz=n*n;return 0;}
00347 static int OneMatDot(void* A, double x[], int nn, int n, double *v){
00348   double dtmp=0.0;
00349   int i,j;
00350   for (i=0;i<n;i++){
00351     for (j=0;j<i;j++,x++){dtmp+= (*x);}
00352     dtmp+= (*x);x++;
00353   }
00354   *v=-2*dtmp;
00355   return 0;
00356 }
00357 static int OneMatVecVec(void* A, double x[], int n, double *v){
00358   double dtmp=0.0;
00359   int i;
00360   for (i=0; i<n; i++){dtmp+=x[i];}
00361   *v=-dtmp*dtmp;
00362   return 0;
00363 }
00364 static int OneMatAddMultiple(void*A, double ddd, double vv[], int nn, int n){
00365   int i,j;
00366   for (i=0;i<n;i++){
00367     for (j=0;j<i;j++,vv++){(*vv)+=-ddd;}
00368     (*vv)+=-ddd; vv++;
00369   }
00370   return 0;
00371 }
00372 static int OneMatAddRowMultiple(void*A, int nrow, double ddd, double row[], int n){
00373   int i;
00374   for (i=0;i<n;i++){row[i] -= ddd;}
00375   row[nrow] -= ddd;
00376   return 0;
00377 }
00378 static int OneMatGetEig(void*A, int neig, double *eig, double v[], int n, int spind[], int *nind){
00379   int i;
00380   if (neig==0){ *eig=-1; for (i=0;i<n;i++){ v[i]=1.0; spind[i]=i;} *nind=n; 
00381   } else {  *eig=0; for (i=0;i<n;i++){ v[i]=0.0; } *nind=0;
00382   }
00383   return 0;
00384 }
00385 static int OneMatRowNnz(void*A, int row, int nz[], int *nnz, int n){
00386   int i;
00387   for (i=0;i<n;i++){ nz[i]++; }
00388   *nnz=n;
00389   return 0;
00390 }
00391 static int OneMatView(void* AA){
00392   printf("Every element of the matrix is the same: -1\n");
00393   return 0;
00394 }
00395 static int OneMatDestroy(void* A){return 0;}
00396 
00397 static const char *mat1name="THETA ALL ELEMENTS EQUAL -ONE";
00398 static int OneMatOpsInitialize(struct  DSDPDataMat_Ops* mat1ops){
00399   int info;
00400   if (mat1ops==NULL) return 0;
00401   info=DSDPDataMatOpsInitialize(mat1ops); DSDPCHKERR(info);
00402   mat1ops->matfactor1=OneMatFactor;
00403   mat1ops->matgetrank=OneMatGetRank;
00404   mat1ops->matgeteig=OneMatGetEig;
00405   mat1ops->matvecvec=OneMatVecVec;
00406   mat1ops->matdot=OneMatDot;
00407   mat1ops->mataddrowmultiple=OneMatAddRowMultiple;
00408   mat1ops->mataddallmultiple=OneMatAddMultiple;
00409   mat1ops->matdestroy=OneMatDestroy;
00410   mat1ops->matview=OneMatView;
00411   mat1ops->matrownz=OneMatRowNnz;
00412   mat1ops->matfnorm2=OneMatFNorm2;
00413   mat1ops->matnnz=OneMatCountNonzeros;
00414   mat1ops->id=18;
00415   mat1ops->matname=mat1name;
00416   return 0;
00417 }
00418 
00419 /* Special Matrix for the Identity */
00420 static int IdentityMatDestroy(void*);
00421 static int IdentityMatView(void*);
00422 static int IdentityMatVecVec(void*, double[], int, double *);
00423 static int IdentityMatDot(void*, double[], int, int, double *);
00424 static int IdentityMatGetRank(void*, int*,int);
00425 static int IdentityMatFactor(void*);
00426 static int IdentityMatGetEig(void*, int, double*, double[], int, int[], int*);
00427 static int IdentityMatAddRowMultiple(void*, int, double, double[], int);
00428 static int IdentityMatAddMultiple(void*, double, double[], int, int);
00429 static int IdentityMatGetRowNnz(void*, int, int[], int*, int);
00430 
00431 static int IdentityMatDestroy(void* AA){return 0;}
00432 static int IdentityMatFNorm2(void* AA, int n, double *v){*v=1.0*n;return 0;}
00433 static int IdentityMatGetRank(void *AA, int*rank, int n){*rank=n;return 0;}
00434 static int IdentityMatFactor(void*A){return 0;}
00435 static int IdentityMatCountNonzeros(void*A, int *nnz, int n){*nnz=n;return 0;}
00436 static int IdentityMatVecVec(void* AA, double x[], int n, double *v){
00437   int i;
00438   *v=0;
00439   for (i=0;i<n;i++){ *v+=x[i]*x[i]; }
00440   return 0;
00441 }
00442 static int IdentityMatDot(void* AA, double x[], int nn, int n, double *v){
00443   int i;
00444   double vv=0;
00445   for (i=0;i<n;i++){ vv+=x[((i+1)*(i+2))/2-1];}
00446   *v = 2*vv;
00447   return 0;
00448 }
00449 static int IdentityMatView(void* AA){
00450   printf("Identity matrix: All Diagonal elements equal 1.0\n");
00451   return 0;
00452 }
00453 static int IdentityMatGetEig(void*AA, int neig, double *eig, double v[], int n, int spind[], int *nind){
00454   if (neig<0 || neig>=n){ *eig=0; return 0;} 
00455   memset((void*)v,0,(n)*sizeof(double)); 
00456   *eig=1.0;  v[neig]=1.0; spind[0]=neig; *nind=1;
00457   return 0;
00458 }
00459 static int IdentityMatGetRowNnz(void*A, int nrow, int nz[], int *nnzz, int n){
00460   if (nrow>=0 && nrow < n){
00461     *nnzz=1; nz[nrow]++;
00462   } else { *nnzz=0;    
00463   }
00464   return 0;
00465 }
00466 static int IdentityMatAddRowMultiple(void*A, int nrow, double dd, double rrow[], int n){
00467   rrow[nrow] += dd;return 0;
00468 }
00469 static int IdentityMatAddMultiple(void*A, double dd, double vv[], int nn, int n){
00470   int i;
00471   for (i=0;i<n;i++){ vv[(i+1)*(i+2)/2-1] += dd;}
00472   return 0;
00473 }
00474 
00475 static const char *eyematname="THETA IDENTITY MATRIX";
00476 static int IdentitymatOperationsInitialize(struct  DSDPDataMat_Ops* imatops){
00477   int info;
00478   if (imatops==NULL) return 0;
00479   info=DSDPDataMatOpsInitialize(imatops); if (info){return info;}
00480   imatops->matfactor1=IdentityMatFactor;
00481   imatops->matgetrank=IdentityMatGetRank;
00482   imatops->matgeteig=IdentityMatGetEig;
00483   imatops->matvecvec=IdentityMatVecVec;
00484   imatops->matrownz=IdentityMatGetRowNnz;
00485   imatops->matdot=IdentityMatDot;
00486   imatops->matfnorm2=IdentityMatFNorm2;
00487   imatops->matnnz=IdentityMatCountNonzeros;
00488   imatops->mataddrowmultiple=IdentityMatAddRowMultiple;
00489   imatops->mataddallmultiple=IdentityMatAddMultiple;
00490   imatops->matdestroy=IdentityMatDestroy;
00491   imatops->matview=IdentityMatView;
00492   imatops->id=12;
00493   imatops->matname=eyematname;
00494   return 0;
00495 }

Generated on Sun Mar 23 07:30:49 2008 for DSDP by  doxygen 1.5.5