DSDP
|
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 }