DSDP
|
00001 #include "dsdp5.h" 00002 #include <string.h> 00003 #include <stdio.h> 00004 #include <math.h> 00005 #include <stdlib.h> 00006 00011 static char help[]="\ 00012 DSDP Usage: dsdp5 filename <sdpa format file> \n\ 00013 -print <10> - print information at each k iteration\n\ 00014 -save <Solution File in SDPA format> \n\ 00015 -fout <filename> to print standard monitor to a file\n\ 00016 -y0 <initial solution file> \n\ 00017 -benchmark <file containing names of SDPA files> \n\ 00018 -directory <directory containing benchmark SDPA files> \n\ 00019 -suffix <add to each benchmark problem name> \n\ 00020 -dloginfo <0> - print more information for higher numbers \n\ 00021 -dlogsummary <1> - print timing information \n\ 00022 -help for this help message\n"; 00023 00024 static int qusort(int[],int[],int[],double[],int,int); 00025 static int partition(int[],int[],int[],double[],int, int); 00026 static int Parseline(char *,int *,int *,int *,int *,double *, int *); 00027 static int ReadInitialPoint(char*, int, double[]); 00028 static int TCheckArgs0(DSDP,SDPCone,int,int,char *[]); 00029 static int TCheckArgs(DSDP,SDPCone,int,int,char *[]); 00030 static int CheckForConstantMat(double[],int, int); 00031 static int CountNonzeroMatrices(int, int[],int[], int*); 00032 00033 typedef struct{ 00034 char sformat; 00035 int blocksize; 00036 } DBlock; 00037 00038 typedef struct{ 00039 int *block,*constraint,*matind; 00040 double*nnz; 00041 char *sformat; 00042 int totalnonzeros; 00043 double *dobj,*y0; 00044 char *conetypes; 00045 int *blocksizes; 00046 int m; int n; int nblocks; 00047 int lpn,lpspot,lpblock,lpnnz; 00048 int *lpi,*lui,*cmap; 00049 double cnorm; 00050 double fixedvari; 00051 double fixedvard,xout; 00052 } DSDPData; 00053 00054 static int ReadSDPA2(char*,DSDPData*); 00055 static int GetMarkers(int, int, int*, int*, int*); 00056 static int ComputeY0(DSDP,DSDPData); 00057 static int rank=0; 00058 int ReadSDPAFile(int argc,char *argv[]); 00059 00060 #define CHKDATA(a,b,c) { if (c){ printf("Possible problem in variable %d, block %d. \n",a+1,b+1);} } 00061 00062 00063 #undef __FUNCT__ 00064 #define __FUNCT__ "main" 00065 int main(int argc,char *argv[]){ 00066 int info; 00067 info=ReadSDPAFile(argc,argv); 00068 return info; 00069 } 00070 00071 #undef __FUNCT__ 00072 #define __FUNCT__ "ReadSDPAFile" 00073 00080 int ReadSDPAFile(int argc,char *argv[]){ 00081 00082 int i,j,m,n,np,its,info; 00083 int spot,ijnnz,nzmats,sdpnmax,sdpn,stat1,printsummary=1; 00084 int runbenchmark=0,saveit=0,justone=1,fileout=0; 00085 double t1,t2,t3,t4,t5,dd,yhigh; 00086 double derr[6],dnorm[3]; 00087 double ddobj,ppobj,scl,dpot; 00088 char problemname[100],thisline[100], filename[300],savefile[100]; 00089 char directory[100]="/home/benson/sdpexamples/sdplib/"; 00090 char outputfile[50]="",suffix[20]=".dat-s", tablename[20]="results-dsdp-5.8"; 00091 char success='s',sformat; 00092 FILE *fp1=0,*fp2=0,*fout; 00093 DSDPData dddd; 00094 DSDP dsdp; 00095 DSDPTerminationReason reason; 00096 DSDPSolutionType pdfeasible; 00097 SDPCone sdpcone=0; 00098 LPCone lpcone=0; 00099 int *ittt,sspot,ione=1; 00100 00101 if (argc<2){ printf("%s",help); DSDPPrintOptions(); return(0); } 00102 for (i=0; i<argc; i++){ if (strncmp(argv[i],"-help",5)==0){ 00103 printf("%s",help); DSDPPrintOptions(); return(0);}} 00104 for (i=1; i<argc; i++){ printf("%s ",argv[i]); } printf("\n"); 00105 for (i=1; i<argc-1; i++){ /* Are we reading a file with a lot of problem names? */ 00106 if (strncmp(argv[i],"-benchmark",8)==0){ 00107 strncpy(thisline,argv[i+1],90); fp1=fopen(thisline,"r");runbenchmark=1; justone=0; 00108 }; 00109 if (strncmp(argv[i],"-directory",8)==0){strncpy(directory,argv[i+1],90);} 00110 if (strncmp(argv[i],"-table",4)==0){strncpy(tablename,argv[i+1],90);}; 00111 if (strncmp(argv[i],"-suffix",4)==0){strncpy(suffix,argv[i+1],20);}; 00112 if (strncmp(argv[i],"-save",5)==0){ strncpy(savefile,argv[i+1],40);saveit=1;}; 00113 if (strncmp(argv[i],"-dlogsummary",8)==0){printsummary=atoi(argv[i+1]);} 00114 if (rank==0&&strncmp(argv[i],"-fout",5)==0){ strncpy(outputfile,argv[i+1],45);fileout=1;}; 00115 } 00116 00117 if (runbenchmark || argc>2){ 00118 fp2=fopen(tablename,"a"); 00119 for (i=1; i<argc; i++){ fprintf(fp2,"%s ",argv[i]); } fprintf(fp2,"\n"); 00120 fclose(fp2); 00121 } 00122 00123 while ((runbenchmark && !feof(fp1)) || justone==1){ 00124 justone=0; 00125 if (runbenchmark){ /* Get filename with the data */ 00126 fgets(thisline,100,fp1); if (sscanf(thisline,"%s",problemname)<1) break; 00127 strncpy(filename,directory,90); strncat(filename,problemname,90);strncat(filename,suffix,90); 00128 printf("%s\n",problemname); 00129 } else { 00130 strncpy(filename,argv[1],90); 00131 strncpy(problemname,argv[1],90); 00132 } 00133 00134 if (rank==0 && fileout){ 00135 dsdpoutputfile=fopen(outputfile,"a"); 00136 fprintf(dsdpoutputfile,"%s\n",problemname); 00137 } else {dsdpoutputfile=0;fileout=0;} 00138 DSDPTime(&t1); 00139 00140 info=ReadSDPA2(filename, &dddd); m=dddd.m; 00141 if (info){ printf("Problem reading SDPA file\n"); return 1;} 00142 DSDPTime(&t2); 00143 if (printsummary && rank==0){ 00144 printf("\nVariables %d \n",dddd.m); 00145 printf("Matrix Blocks: %d, ",dddd.nblocks); 00146 printf("Total Number of Constraints: %d \n",dddd.n); 00147 printf("Nonzeros in Constraints: %d\n\n",dddd.totalnonzeros); 00148 printf("Read Data File into Buffer: %4.3e seconds\n",t2-t1); 00149 } 00150 00151 for (i=0; i<argc-1; i++){ 00152 if (strncmp(argv[i],"-dloginfo",8)==0){ info=DSDPLogInfoAllow(atoi(argv[i+1]),0);}; 00153 } 00154 00155 info = DSDPCreate(dddd.m,&dsdp); 00156 info = DSDPCreateSDPCone(dsdp,dddd.nblocks,&sdpcone); 00157 /* Set Dual objective vector */ 00158 for (i=0;i<m;i++){info = DSDPSetDualObjective(dsdp,i+1,dddd.dobj[i]);} 00159 00160 /* Set initial point */ 00161 for (i=0; i<m; i++) 00162 if (dddd.dobj[i]> 0.0) dddd.y0[i]=-0.0; else dddd.y0[i]=0.0; 00163 for (i=0; i<m; i++){info = DSDPSetY0(dsdp,i+1,dddd.y0[i]);} 00164 info=ComputeY0(dsdp,dddd); 00165 if (dddd.fixedvari){ 00166 info = DSDPSetY0(dsdp,(int)dddd.fixedvari,dddd.fixedvard); 00167 printf("Fixed: %2.0f %4.2f ?\n",dddd.fixedvari,dddd.fixedvard); 00168 info=DSDPSetFixedVariables(dsdp,&dddd.fixedvari,&dddd.fixedvard,&dddd.xout,ione); 00169 } 00170 00171 spot=0;ijnnz=0;np=0;sdpnmax=1;sdpn=0;stat1=1; 00172 /* Insert the SDP data */ 00173 for (j=0;j<dddd.nblocks; j++){ 00174 if (dddd.conetypes[j]=='S'){ 00175 n=dddd.blocksizes[j]; 00176 sformat=dddd.sformat[j]; 00177 info=CountNonzeroMatrices(j+1,dddd.block+spot,dddd.constraint+spot,&nzmats); 00178 info=SDPConeSetBlockSize(sdpcone,j,n); DSDPCHKERR(info); 00179 info=SDPConeSetSparsity(sdpcone,j,nzmats); DSDPCHKERR(info); 00180 info=SDPConeSetStorageFormat(sdpcone,j,sformat); DSDPCHKERR(info); 00181 np+=n; sdpn+=n; 00182 if (sdpnmax<n) sdpnmax=n; 00183 if (stat1<nzmats) stat1=nzmats; 00184 for (i=0; i<=m; i++){ 00185 info=GetMarkers(j+1,i,dddd.block+spot,dddd.constraint+spot,&ijnnz); 00186 if (0==1){ 00187 } else if ( ijnnz==0 ){ /* info=DSDPSetZeroMat(dsdp,j,i,n); */ 00188 } else if (CheckForConstantMat(dddd.nnz+spot,ijnnz,n)){ 00189 info=SDPConeSetConstantMat(sdpcone,j,i,n,dddd.nnz[spot+1]);CHKDATA(i,j,info); 00190 if(sformat=='P'){info=SDPConeSetXArray(sdpcone,j,n,dddd.nnz+spot,n*(n+1)/2);} 00191 } else if (sformat=='P' && ijnnz==n*(n+1)/2 ){ /* check for dense matrix */ 00192 info=SDPConeSetADenseVecMat(sdpcone,j,i,n,1.0,dddd.nnz+spot,ijnnz);CHKDATA(i,j,info); 00193 } else { /* sparse matrix */ 00194 info=SDPConeSetASparseVecMat(sdpcone,j,i,n,1.0,0,dddd.matind+spot,dddd.nnz+spot,ijnnz);CHKDATA(i,j,info); 00195 } 00196 if (0==1){ info=SDPConeViewDataMatrix(sdpcone,j,i);} 00197 spot+=ijnnz; 00198 /* SDPConeScaleBarrier(sdpcone,j,j+1.0); */ 00199 } 00200 if (0==1){info=SDPConeView2(sdpcone);} 00201 } else if (dddd.conetypes[j]=='L'){ 00202 info=DSDPCreateLPCone(dsdp,&lpcone); sformat='P'; 00203 info=SDPConeSetStorageFormat(sdpcone,j,sformat); DSDPCHKERR(info); 00204 n=dddd.blocksizes[j]; 00205 np+=n; 00206 sspot=spot; 00207 DSDPCALLOC2(&ittt,int,(m+2),&info); 00208 for (i=0;i<=m;i++){ittt[i]=0;} 00209 for (i=0;i<=m;i++){ 00210 info=GetMarkers(j+1,i,dddd.block+spot,dddd.constraint+spot,&ijnnz); 00211 ittt[i+1]=ijnnz; spot+=ijnnz; 00212 } 00213 for (i=1;i<=m;i++)ittt[i+1]+=ittt[i]; 00214 info=LPConeSetData(lpcone,n,ittt,dddd.matind+sspot,dddd.nnz+sspot);CHKDATA(i,0,info); 00215 if (0==1){info=LPConeView(lpcone);} 00216 if (0==1){info=LPConeView2(lpcone);} 00217 /* info=DSDPFree(&ittt); */ 00218 } 00219 } 00220 if (0==1){ 00221 BCone bcone; 00222 info=DSDPCreateBCone(dsdp, &bcone); 00223 info=BConeAllocateBounds(bcone,2*m); 00224 for (i=0;i<m;i++){ 00225 info=BConeSetUpperBound(bcone,i+1,10); 00226 } 00227 for (i=0;i<m;i++){ 00228 info=BConeSetLowerBound(bcone,i+1,-10); 00229 } 00230 } 00231 00232 DSDPTime(&t3); 00233 if (printsummary && rank==0){printf("DSDP Set Data: %4.3e seconds\n",t3-t2);} 00234 00235 its=(m-2)/sdpnmax; 00236 if (np<100 && its==0) its=1; 00237 if (its>=1) its++; 00238 its=its*its; 00239 if (m<2000 && its>10) its=10; 00240 if (its>12) its=12; 00241 00242 info=DSDPReuseMatrix(dsdp,its); 00243 00244 DSDPFREE(&dddd.blocksizes,&info); 00245 DSDPFREE(&dddd.sformat,&info); 00246 DSDPFREE(&dddd.dobj,&info); 00247 DSDPFREE(&dddd.y0,&info); 00248 DSDPFREE(&dddd.conetypes,&info); 00249 DSDPFREE(&dddd.constraint,&info); 00250 DSDPFREE(&dddd.block,&info); 00251 00252 info=DSDPGetDataNorms(dsdp, dnorm); 00253 if (dnorm[0]==0){ 00254 info=DSDPSetR0(dsdp,np); 00255 info=DSDPSetGapTolerance(dsdp,1e-3); 00256 info=DSDPSetYBounds(dsdp,-1.0,1.0); 00257 } else { 00258 } 00259 info = TCheckArgs0(dsdp,sdpcone,dddd.m,argc,argv); 00260 info = TCheckArgs(dsdp,sdpcone,dddd.m,argc,argv); 00261 00262 info = DSDPSetup(dsdp); if (info){ printf("\nProblem Setting problem. Likely insufficient memory\n"); return 1;} 00263 if (0==1){info=SDPConeCheckData(sdpcone);} 00264 00265 00266 DSDPTime(&t4); 00267 info=DSDPGetScale(dsdp,&scl); 00268 info=DSDPGetPotentialParameter(dsdp,&dpot); 00269 info=DSDPGetReuseMatrix(dsdp,&its); 00270 if (printsummary && rank==0){ 00271 printf("DSDP Process Data: %4.3e seconds\n\n",t4-t3); 00272 printf("Data Norms: C: %4.2e, A: %4.2e, b: %4.2e\n",dnorm[0],dnorm[1],dnorm[2]); 00273 printf("Scale C: %4.2e\n\n",scl); 00274 printf("Potential Parameter: %4.2f\n",dpot); 00275 printf("Reapply Schur matrix: %d\n\n",its); 00276 } 00277 if (0==1){info=DSDPPrintData(dsdp,sdpcone,lpcone);} 00278 00279 info = DSDPSolve(dsdp); 00280 if (info){ printf("\nNumerical errors encountered in DSDPSolve(). \n");} 00281 00282 info=DSDPStopReason(dsdp,&reason); 00283 if (reason!=DSDP_INFEASIBLE_START){ 00284 info=DSDPComputeX(dsdp);DSDPCHKERR(info); 00285 } 00286 info=DSDPStopReason(dsdp,&reason); 00287 info=DSDPGetSolutionType(dsdp,&pdfeasible); 00288 00289 DSDPTime(&t5); 00290 00291 info=DSDPGetDObjective(dsdp,&ddobj); 00292 info=DSDPGetPObjective(dsdp,&ppobj); 00293 info=DSDPGetFinalErrors(dsdp,derr); 00294 info=DSDPGetIts(dsdp,&its); 00295 00296 success='s'; 00297 if (printsummary && rank==0){ 00298 00299 if (reason == DSDP_CONVERGED){ 00300 printf("DSDP Converged. \n"); 00301 success='s'; 00302 } else if ( reason == DSDP_UPPERBOUND ){ 00303 printf("DSDP Terminated Because Dual Objective Exceeded its Bound\n"); 00304 success='c'; 00305 } else if ( reason == DSDP_SMALL_STEPS ){ 00306 printf("DSDP Terminated Due to Small Steps\n"); 00307 success='c'; 00308 } else if ( reason == DSDP_MAX_IT){ 00309 printf("DSDP Terminated Due Maximum Number of Iterations\n"); 00310 success='c'; 00311 } else if ( reason == DSDP_INFEASIBLE_START){ 00312 printf("DSDP Terminated Due to Infeasible Starting Point\n"); 00313 success='c'; 00314 } else if ( reason == DSDP_INDEFINITE_SCHUR_MATRIX){ 00315 printf("DSDP Terminated Due to Indefinite Schur Complement\n"); 00316 success='c'; 00317 } else { 00318 printf("DSDP Finished\n"); 00319 success='c'; 00320 } 00321 00322 if (pdfeasible == DSDP_UNBOUNDED ){ 00323 printf("DSDP Dual Unbounded, Primal Infeasible\n"); 00324 } else if ( pdfeasible == DSDP_INFEASIBLE ){ 00325 printf("DSDP Primal Unbounded, Dual Infeasible\n"); 00326 } 00327 00328 00329 printf("\nP Objective : %16.8e \n",ppobj); 00330 printf("DSDP Solution: %16.8e \n\n",ddobj); 00331 printf("DSDP Solve Time: %4.3e seconds\n",t5-t4); 00332 printf("DSDP Preparation and Solve Time: %4.3e seconds\n\n",t5-t3); 00333 00334 } else { 00335 if (reason == DSDP_CONVERGED){success='s';} else {success='c';} 00336 } 00337 00338 if (rank==0){ 00339 fp2=fopen(tablename,"a"); 00340 if (pdfeasible==DSDP_UNBOUNDED){ 00341 fprintf(fp2," %-18s & %4d & %4d & infeasible & unbounded & %4.0e & %c & %3d & %6.2f \\\\\n",problemname,m,np,derr[0],success,its,t5-t3); 00342 } else if (pdfeasible==DSDP_INFEASIBLE){ 00343 fprintf(fp2," %-18s & %4d & %4d & unbounded & infeasible & %4.0e & %c & %3d & %6.2f \\\\\n",problemname,m,np,derr[0],success,its,t5-t3); 00344 } else { 00345 fprintf(fp2," %-18s & %4d & %4d & %13.7f & %13.7f & %4.0e & %c & %3d & %6.2f \\\\\n",problemname,m,np,-ppobj,-ddobj,derr[0],success,its,t5-t3); 00346 } 00347 fclose(fp2); 00348 } 00349 00350 if (printsummary && rank==0){ 00351 /* info=DSDPComputeMinimumXEigenvalue(dsdp,&derr[1]); */ 00352 printf("\nP Infeasible: %8.2e \n",derr[0]); 00353 printf("D Infeasible: %8.2e \n",derr[2]); 00354 printf("Minimal P Eigenvalue: %6.2e \n",derr[1]); 00355 printf("Minimal D Eigenvalue: 0.00 \n"); 00356 printf("Relative P - D Objective values: %4.2e \n",derr[4]); 00357 printf("Relative X Dot S: %4.2e \n",derr[5]); 00358 info=DSDPGetYBounds(dsdp,&dd,&yhigh); 00359 info=DSDPGetYMaxNorm(dsdp,&dd); 00360 printf("\nMax Y: %10.8e, Bounded by %6.1e\n",dd,yhigh); 00361 info=DSDPGetTraceX(dsdp,&dd); 00362 printf("Trace X: %4.8e, ",dd); 00363 info=DSDPGetPenaltyParameter(dsdp,&dd); 00364 printf("Bounded by Penalty Parameter: %4.1e \n\n",dd); 00365 00366 if (printsummary){ DSDPEventLogSummary();} 00367 printf("--- DSDP Finished ---\n\n"); 00368 } 00369 00370 if (rank==0){ 00371 if (saveit){ 00372 fout=fopen(savefile,"w"); 00373 /* fprintf(fout,"** %s \n",filename); Deleted */ 00374 info= DSDPPrintSolution(fout,dsdp,sdpcone,lpcone); 00375 if (dddd.fixedvari){ 00376 sspot=dddd.nblocks+1,dd=dddd.xout; 00377 fprintf(fout,"1 %d 1 1 1.0e-11\n1 %d 2 2 1.0e-11\n",sspot,sspot); 00378 fprintf(fout,"2 %d 1 1 %12.8e\n",sspot,DSDPMax(1.0e-10,dd)); 00379 fprintf(fout,"2 %d 2 2 %12.8e\n",sspot,DSDPMax(1e-10,-dd)); 00380 } 00381 fclose(fout); 00382 } 00383 } 00384 00385 for (i=0; i<argc; i++){ 00386 if (rank==0 && strncmp(argv[i],"-view",5)==0){DSDPView(dsdp);} 00387 } 00388 00389 info = DSDPDestroy(dsdp); 00390 DSDPFREE(&dddd.matind,&info); 00391 DSDPFREE(&dddd.nnz,&info); 00392 00393 if (fileout){fclose(dsdpoutputfile);} 00394 if (0){ DSDPMemoryLog();} 00395 00396 } 00397 if (runbenchmark){ fclose(fp1);} 00398 00399 return 0; 00400 } /* main */ 00401 00402 00403 00404 #define BUFFERSIZ 4000 00405 #define BUFFERSIZ 4000 00406 #undef __FUNCT__ 00407 #define __FUNCT__ "ReadSDPA2" 00408 static int ReadSDPA2(char *filename, DSDPData*ddd){ 00409 FILE*fp; 00410 char ctmp,refline[BUFFERSIZ]="*",thisline[BUFFERSIZ]="*"; 00411 int info,tline,line=0; 00412 int i,k,m,n; 00413 /* int spot,nzmark, */ 00414 int bigint=1000000; 00415 int i1,nblk,nmat,col,row; 00416 int np=0,nblocks; 00417 int nargs,nonzero; 00418 double val; 00419 00420 fp=fopen(filename,"r"); 00421 if (!fp){ 00422 printf("Cannot open file %s !",filename); return(1); 00423 } 00424 00425 /* Read comments */ 00426 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){ 00427 fgets(thisline,BUFFERSIZ,fp); line++; 00428 } 00429 /* Read number of constraints */ 00430 if (sscanf(thisline,"%d",&m)<1){ 00431 printf("Error: line %d. Number of constraints not given.\n",line); 00432 return(1); 00433 } 00434 /* Read number of blocks */ 00435 fgets(thisline,BUFFERSIZ,fp); line++; 00436 if (sscanf(thisline,"%d",&nblocks)!=1){ 00437 printf("Error: line %d. Number of blocks not given.\n",line); 00438 return(1); 00439 } 00440 ddd->lpn=0;ddd->lpspot=0;ddd->lpblock=0;ddd->cnorm=0; 00441 /* Read block sizes */ 00442 DSDPCALLOC2(&ddd->sformat,char, (nblocks+1),&info); 00443 DSDPCALLOC2(&ddd->blocksizes,int, (nblocks+1),&info); 00444 DSDPCALLOC2(&ddd->conetypes,char, (nblocks+1),&info ); 00445 line++; 00446 for (i=0;i<nblocks; i++){ 00447 if (fscanf(fp,"{")==1 || fscanf(fp,"(")==1 || fscanf(fp,",")==1 ){ 00448 i--; 00449 } else if (fscanf(fp,"%d",&col)==1){ 00450 if (col>0) { ddd->blocksizes[i]=col; np+=col; ddd->conetypes[i]='S'; 00451 } else if (col>0){ddd->blocksizes[i]=-col; np+=-col; ddd->conetypes[i]='S'; 00452 } else if (col<0){ddd->blocksizes[i]=-col; np += -col; ddd->conetypes[i]='L';ddd->lpn=-col;ddd->lpblock=i; 00453 } else { ddd->blocksizes[i]=0; ddd->conetypes[i]='N';} 00454 if (ddd->blocksizes[i]<10){ddd->sformat[i]='U';} else {ddd->sformat[i]='P';} 00455 } 00456 else{ printf("Error block sizes, line %d",line); return(1);} 00457 } 00458 if (ddd->blocksizes[nblocks-1]==0) nblocks--; 00459 fgets(thisline,BUFFERSIZ,fp); 00460 00461 /* Read objective vector */ 00462 DSDPCALLOC2(&ddd->y0,double,m,&info); 00463 DSDPCALLOC2(&ddd->dobj,double,m,&info); 00464 line++; 00465 for (i=0;i<m;i++){ 00466 if (fscanf(fp,",")==1){ 00467 i--; 00468 continue; 00469 } 00470 while (fscanf(fp,"%lg",&val)!=1){ 00471 fscanf(fp,"%c",&ctmp); 00472 if (ctmp=='\n'){ 00473 printf("Constraints: %d, Blocks: %d\n",m,nblocks); 00474 printf("Error reading objective, line %d, i=%d \n",line,i); return 1; 00475 } 00476 } 00477 ddd->dobj[i]=val; 00478 } 00479 00480 fgets(thisline,BUFFERSIZ,fp); 00481 tline=line; 00482 00483 fseek(fp,0,SEEK_SET); 00484 line=0; 00485 for (i=0;i<tline;i++){ctmp='*'; while (ctmp!='\n'){ fscanf(fp,"%c",&ctmp);} line++;} 00486 00487 nargs=5; nonzero=0; 00488 while(!feof(fp)){ 00489 thisline[0]='\0'; 00490 nmat=-1; nblk=-1; row=-1; col=-1; val=0.0; 00491 fgets(thisline,BUFFERSIZ,fp); line++; 00492 info = Parseline(thisline,&nmat,&nblk,&row,&col,&val,&nargs); 00493 if (!feof(fp)&&nargs!=5&&nargs>0){ 00494 printf("Error: line %d \n%s\n",line,thisline);return 1;} 00495 if (nargs==5 && val!=0.0){ 00496 nonzero++; 00497 i1=row*(row+1)/2 + col; 00498 if (row >= ddd->blocksizes[nblk-1] || col >= ddd->blocksizes[nblk-1] ) { 00499 printf("Data Error in line: %d. Row %d or col %d > blocksize %d\n%s",line,row+1,col+1,ddd->blocksizes[nblk-1],thisline); 00500 return 1; 00501 } 00502 if (row<0 || col<0){ 00503 printf("Data Error in line: %d. Row %d or col %d <= 0 \n%s",line,row+1,col+1,thisline); 00504 return 1; 00505 } 00506 if (nmat>m || nmat<0){ 00507 printf("Data Error in line: %d. Is Var 0 <= %d <= %d \n%s",line,nmat,m,thisline); 00508 return 1; 00509 } 00510 if (nblk>nblocks || nblk<0){ 00511 printf("Data Error in line: %d. Is Block 0 <= %d <= %d \n%s",line,nmat,m,thisline); 00512 return 1; 00513 } 00514 } else if (nargs==5 && val==0.0){ 00515 } else if (nargs<5 && nargs>0){ 00516 printf("Too few numbers. \n"); 00517 printf("Problem Reading SDPA file at line %d: %s\n",line, thisline); 00518 } else if (nargs==0){ 00519 } else { 00520 printf("Problem Reading SDPA file at line %d: %s\n",line, thisline); 00521 } 00522 } 00523 00524 /* Allocate memory for the data */ 00525 nonzero++; 00526 DSDPCALLOC2(&ddd->matind,int,nonzero,&info); 00527 DSDPCALLOC2(&ddd->nnz,double,nonzero,&info); 00528 DSDPCALLOC2(&ddd->block,int,nonzero,&info); 00529 DSDPCALLOC2(&ddd->constraint,int,nonzero,&info); 00530 nonzero--; 00531 00532 fseek(fp,0,SEEK_SET); 00533 line=0; 00534 for (i=0;i<tline;i++){ctmp='*'; while (ctmp!='\n') fscanf(fp,"%c",&ctmp); line++;} 00535 00536 nargs=5;k=0; 00537 while(!feof(fp) /* && nargs==5 */){ 00538 thisline[0]='\0'; 00539 fgets(thisline,BUFFERSIZ,fp); 00540 if (k==0){strncpy(refline,thisline,BUFFERSIZ-1); } 00541 info = Parseline(thisline,&nmat,&nblk,&row,&col,&val,&nargs); 00542 if (!feof(fp)&&nargs!=5&&nargs<0){ 00543 /* if (k>=nonzero && !feof(fp) ){ */ 00544 printf("Problem Reading SDPA file at line %d: %s\n",line, thisline); 00545 printf("Problem could be earlier in file \n"); 00546 printf("The first recorded matix nonzero in the file occured at line %d: \n %s",tline,refline); 00547 printf(" Please check data file\n"); 00548 return 1; 00549 } 00550 if (nargs==5 && val!=0.0){ 00551 if (row>col){ 00552 printf("Warning: Line: %d Row < Column. %s \n",line,thisline); 00553 } 00554 i=row;row=col;col=i; 00555 n=ddd->blocksizes[nblk-1]; 00556 if (nmat==0) {val=-val;} 00557 if (ddd->conetypes[nblk-1]=='S'){ 00558 /* if (row==col) val/=2; */ 00559 ddd->matind[k]=row*(row+1)/2 + col; 00560 if (ddd->sformat[nblk-1]=='U'){ddd->matind[k]=row*n + col;} 00561 } else { 00562 ddd->matind[k]=col; 00563 } 00564 ddd->block[k]=nblk; 00565 ddd->constraint[k]=nmat; 00566 ddd->nnz[k]=val; 00567 k++; 00568 } else if (nargs==5 && val==0.0){ 00569 } else if (nargs==0){ 00570 } else { 00571 printf("Problem Reading SDPA file at line %d: %s\n",line, thisline); 00572 printf("Problem could be earlier in file \n"); 00573 printf("The first recorded matix nonzero in the file occured at line %d: \n %s",tline,refline); 00574 printf(" Please check data file\n"); 00575 return 1; 00576 } 00577 } 00578 ddd->block[k]=nblocks+1; ddd->constraint[k]=m+2; 00579 ddd->matind[k]=10000000; ddd->nnz[k]=0.0; 00580 00581 qusort(ddd->block,ddd->constraint,ddd->matind,ddd->nnz,0,nonzero-1); 00582 00583 for (i=0;i<nonzero-1; i++){ 00584 while (i<nonzero-1 && ddd->matind[i]==ddd->matind[i+1] && 00585 ddd->constraint[i]==ddd->constraint[i+1] && 00586 ddd->block[i]==ddd->block[i+1] ){ 00587 printf("DSDPError: Reading Input File:\n"); 00588 printf("Possible problem with data input file: Double Entry: \n"); 00589 printf(" %d %d %d %2.10e\n", 00590 ddd->constraint[i],ddd->block[i],ddd->matind[i]+1,ddd->nnz[i]); 00591 printf(" %d %d %d %2.10e\n\n", 00592 ddd->constraint[i+1],ddd->block[i+1],ddd->matind[i+1]+1,ddd->nnz[i+1]); 00593 for (k=i+1;k<nonzero-1;k++){ 00594 ddd->constraint[k]=ddd->constraint[k+1]; ddd->block[k]=ddd->block[k+1]; 00595 ddd->matind[k]=ddd->matind[k+1];ddd->nnz[k]=ddd->nnz[k+1]; 00596 } 00597 ddd->constraint[nonzero-1]=bigint;ddd->nnz[nonzero-1]=0; 00598 nonzero--; 00599 } 00600 } 00601 00602 ddd->fixedvari=0;ddd->fixedvard=0; 00603 if (ddd->lpblock>0){ 00604 int spot; 00605 if (ddd->blocksizes[ddd->lpblock]==2){ 00606 i=0;k=0; 00607 while (ddd->block[i]<=ddd->lpblock && i<nonzero){ i++;} spot=i; 00608 while (ddd->block[i]==ddd->lpblock+1 && i<nonzero){ i++;k++;} 00609 if (k==4){ 00610 if (ddd->constraint[spot]==ddd->constraint[spot+1] && 00611 ddd->constraint[spot+2]==ddd->constraint[spot+3] && 00612 ddd->matind[spot]==ddd->matind[spot+2] && 00613 ddd->matind[spot+1]==ddd->matind[spot+3] && 00614 fabs(ddd->nnz[spot+2])==1.0 && fabs(ddd->nnz[spot+3])==1.0 && 00615 fabs(ddd->nnz[spot] + ddd->nnz[spot+1]) <=1e-6 ){ 00616 ddd->fixedvari=ddd->constraint[spot+2]; 00617 ddd->fixedvard=ddd->nnz[spot]/ddd->nnz[spot+2]; 00618 nblocks--;ddd->lpblock=0; 00619 } 00620 } 00621 } 00622 } 00623 00624 ddd->totalnonzeros=nonzero; 00625 for (ddd->n=0,i=0;i<nblocks;i++) ddd->n += ddd->blocksizes[i]; 00626 ddd->m=m; 00627 ddd->nblocks=nblocks; 00628 fclose(fp); 00629 return 0; 00630 } 00631 00632 00633 #undef __FUNCT__ 00634 #define __FUNCT__ "Parseline" 00635 static int Parseline(char thisline[],int *nmat,int *nblk,int *row, 00636 int *col,double *value, int *nargs){ 00637 00638 int temp; 00639 int rtmp,coltmp; 00640 00641 *nargs=0; 00642 *nmat=-1;*nblk=-1;rtmp=-1;coltmp=-1;*value=0.0; 00643 temp=sscanf(thisline,"%d %d %d %d %lg",nmat,nblk,&rtmp,&coltmp,value); 00644 if (temp==5) *nargs=5; 00645 else if (temp>0 && temp<5) *nargs=temp; 00646 *row=rtmp-1; *col=coltmp-1; 00647 00648 return(0); 00649 } 00650 00651 00652 static int partition(int list1[], int list2[], int list3[], double list5[], int lstart, int lend){ 00653 int k=lend; 00654 int pivot1=list1[k],pivot2=list2[k],pivot3=list3[k]; 00655 double pivot5 = list5[k]; 00656 int bottom = lstart-1, top = lend; 00657 int done = 0; 00658 int ordered=1; 00659 while (!done){ 00660 00661 while (!done) { 00662 bottom = bottom+1; 00663 00664 if (bottom == top){ 00665 done = 1; 00666 break; 00667 } 00668 if ( list1[bottom] > pivot1 || 00669 (list1[bottom] == pivot1 && list2[bottom] > pivot2) || 00670 (list1[bottom] == pivot1 && list2[bottom] == pivot2 && 00671 list3[bottom] > pivot3) ){ 00672 list1[top] = list1[bottom]; 00673 list2[top] = list2[bottom]; 00674 list3[top] = list3[bottom]; 00675 list5[top] = list5[bottom]; 00676 ordered=0; 00677 break; 00678 } 00679 } 00680 while (!done){ 00681 top = top-1; 00682 00683 if (top == bottom){ 00684 done = 1; 00685 break; 00686 } 00687 if ( list1[top] < pivot1 || 00688 (list1[top] == pivot1 && list2[top] < pivot2) || 00689 (list1[top] == pivot1 && list2[top] == pivot2 && list3[top] < pivot3)){ 00690 list1[bottom] = list1[top]; 00691 list2[bottom] = list2[top]; 00692 list3[bottom] = list3[top]; 00693 list5[bottom] = list5[top]; 00694 ordered=0; 00695 break; 00696 } 00697 } 00698 } 00699 list1[top] = pivot1; 00700 list2[top] = pivot2; 00701 list3[top] = pivot3; 00702 list5[top] = pivot5; 00703 00704 ordered=0; 00705 00706 if (bottom==lend){ 00707 ordered=1; 00708 for (k=lstart;k<lend-1;k++){ 00709 if ( (list1[k] > list1[k+1]) || 00710 (list1[k] == list1[k+1] && list2[k] > list2[k+1]) || 00711 (list1[k] == list1[k+1] && list2[k] == list2[k+1] && list3[k] > list3[k+1]) ){ 00712 ordered=0; 00713 break; 00714 } 00715 } 00716 } 00717 if (ordered && lend-lstart>2){ 00718 top=(lend+lstart)/2; 00719 } 00720 return top; 00721 } 00722 00723 00724 00725 static int qusort(int list1[], int list2[], int list3[], double list5[], int lstart, int lend){ 00726 int split; 00727 if (lstart < lend){ 00728 split = partition(list1, list2, list3, list5,lstart, lend); 00729 qusort(list1, list2, list3, list5, lstart, split-1); 00730 qusort(list1, list2, list3, list5, split+1, lend); 00731 } else { 00732 return 0; 00733 } 00734 return 0; 00735 } 00736 00737 00738 00739 #undef __FUNCT__ 00740 #define __FUNCT__ "GetMarkers" 00741 static int GetMarkers(int block, int constraint, int blockn[], int constraintn[], 00742 int*m3){ 00743 int i=0; 00744 while (blockn[i]==block && constraintn[i]==constraint){ i++;} 00745 *m3=i; 00746 return 0; 00747 } 00748 00749 #undef __FUNCT__ 00750 #define __FUNCT__ "CountNonzeroMatrices" 00751 static int CountNonzeroMatrices(int block, int blockn[], int constraintn[], int*m3){ 00752 int i=0,cvar=-1,nnzmats=0; 00753 while (blockn[i]==block){ 00754 if (constraintn[i]>cvar){ cvar=constraintn[i];nnzmats++;} 00755 i++; 00756 } 00757 *m3=nnzmats; 00758 return 0; 00759 } 00760 00761 #undef __FUNCT__ 00762 #define __FUNCT__ "CheckForConstantMat" 00763 static int CheckForConstantMat(double v[],int nnz, int n){ 00764 int i; double vv; 00765 if (n<=1){ return 0; } 00766 if (nnz!=(n*n+n)/2){ return 0; } 00767 for (vv=v[0],i=1;i<nnz;i++){ 00768 if (v[i]!=vv){ return 0;} 00769 } 00770 return 1; 00771 } 00772 00773 static int ComputeY0(DSDP dsdp,DSDPData dddd){ 00774 int i,ii,info,ijnnz=0,spot=0,ddiag=0,diag=0,n=dddd.blocksizes[0],m=dddd.m; 00775 double aa,bb=0,ddmax=0,dd=0,cnorm=0; 00776 char sformat=dddd.sformat[0]; 00777 if (dddd.nblocks>1) return 0; 00778 if (dddd.fixedvari) return 0; 00779 00780 info=GetMarkers(1,0,dddd.block+spot,dddd.constraint+spot,&ijnnz); 00781 for (i=0;i<ijnnz;i++){if (cnorm<fabs(dddd.nnz[i])) cnorm=fabs(dddd.nnz[i]);} 00782 spot+=ijnnz; 00783 00784 for (i=1;i<=m;i++,spot+=ijnnz){ 00785 info=GetMarkers(1,i,dddd.block+spot,dddd.constraint+spot,&ijnnz); 00786 ddiag=0; 00787 if (ijnnz==n){ 00788 dd=dddd.nnz[spot]; ddiag=1;bb=dddd.dobj[i-1]; 00789 for (ii=0;ii<n;ii++){ 00790 if (sformat=='U'){ 00791 if (dddd.matind[spot+ii] != ii*n+ii || dddd.nnz[spot+i]!=dd){ ddiag=0;} 00792 } else if (sformat=='P'){ 00793 if (dddd.matind[spot+ii] != ii*(ii+1)/2+ii || dddd.nnz[spot+i]!=dd){ ddiag=0;} 00794 } 00795 } 00796 } 00797 if (ddiag){ 00798 info=DSDPSetY0(dsdp,i,-100*n*cnorm*dddd.dobj[i-1]); 00799 info=DSDPSetR0(dsdp,0); 00800 info=DSDPSetZBar(dsdp,100*dd/bb*cnorm); 00801 info=DSDPSetPotentialParameter(dsdp,5.0); 00802 } 00803 } 00804 00805 if (m==n){ 00806 spot=0; info=GetMarkers(1,0,dddd.block+spot,dddd.constraint+spot,&ijnnz); spot+=ijnnz; 00807 dd=dddd.nnz[spot]; bb=dddd.dobj[0]; 00808 for (diag=1,i=0;i<n;i++,spot+=ijnnz){ 00809 info=GetMarkers(1,i+1,dddd.block+spot,dddd.constraint+spot,&ijnnz); 00810 dd=dddd.nnz[spot];bb=dddd.dobj[i]; 00811 if (ddmax<bb/dd) ddmax=bb/dd; 00812 if (sformat=='P'){ 00813 if (ijnnz!=1 || bb!=dddd.dobj[i] || dd!=dddd.nnz[spot] || dddd.matind[spot]!= ((i+1)*(i+2))/2-1) { diag=0; } 00814 } else if (sformat=='U'){ 00815 if (ijnnz!=1 || bb!=dddd.dobj[i] || dd!=dddd.nnz[spot] || dddd.matind[spot]!= ((i)*n)+i) { diag=0; } 00816 } 00817 } 00818 if (diag && cnorm*sqrt(1.0*n)<1e6){ 00819 for (i=0;i<n;i++){info = DSDPSetY0(dsdp,i+1,-10*sqrt(1.0*n)*cnorm);} 00820 info=DSDPSetR0(dsdp,0); info=DSDPSetZBar(dsdp,100*ddmax*n*cnorm); info=DSDPSetPotentialParameter(dsdp,5.0); info=DSDPReuseMatrix(dsdp,0); 00821 } 00822 } 00823 if (m==n+1){ 00824 spot=0; info=GetMarkers(1,0,dddd.block+spot,dddd.constraint+spot,&ijnnz); spot+=ijnnz; 00825 dd=dddd.nnz[spot]; bb=dddd.dobj[0];diag=1; 00826 info=GetMarkers(1,1,dddd.block+spot,dddd.constraint+spot,&ijnnz); 00827 if (CheckForConstantMat(dddd.nnz+spot,ijnnz,n)){aa=dddd.nnz[0]; spot+=ijnnz;ii=2;} else {ii=1;} 00828 for (i=0;i<n;i++,spot+=ijnnz){ 00829 info=GetMarkers(1,i+ii,dddd.block+spot,dddd.constraint+spot,&ijnnz); 00830 dd=dddd.nnz[spot];bb=dddd.dobj[i+ii-1]; 00831 if (ddmax<bb/dd) ddmax=bb/dd; 00832 if (sformat=='U'){ 00833 if (ijnnz!=1 || bb!=dddd.dobj[ii-1] || dd!=dddd.nnz[spot] || dddd.matind[spot]!= ((i)*(n))+i ) { diag=0; } 00834 } else if (sformat=='P'){ 00835 if (ijnnz!=1 || bb!=dddd.dobj[ii-1] || dd!=dddd.nnz[spot] || dddd.matind[spot]!= ((i+1)*(i+2))/2-1) { diag=0; } 00836 } 00837 } 00838 if (ii==1 && diag==1){ 00839 info=GetMarkers(1,m,dddd.block+spot,dddd.constraint+spot,&ijnnz); 00840 if (CheckForConstantMat(dddd.nnz+spot,ijnnz,n)){aa=dddd.nnz[spot];} else {diag=0;} 00841 } 00842 if (diag && cnorm*sqrt(1.0*n)<1e5){ 00843 /* 00844 if (ii=2){info = DSDPSetY0(dsdp,1,-10000*aa);} else {info = DSDPSetY0(dsdp,m,-10000*aa);} 00845 for (i=0;i<n;i++){info = DSDPSetY0(dsdp,i+ii,-100*sqrt(1.0*n)*cnorm);} 00846 */ 00847 /* info=DSDPSetR0(dsdp,cnorm); info=DSDPSetZBar(dsdp,n*n*n*ddmax*cnorm); */ info=DSDPSetPotentialParameter(dsdp,5.0); info=DSDPReuseMatrix(dsdp,0); 00848 } 00849 } 00850 return 0; 00851 } 00852 00853 #undef __FUNCT__ 00854 #define __FUNCT__ "TCheckArgs0" 00855 static int TCheckArgs0(DSDP dsdp,SDPCone sdpcone, int m,int nargs,char *runargs[]){ 00856 00857 int kk,info,iloginfo=0; 00858 for (kk=1; kk<nargs-1; kk++){ 00859 if (0){ 00860 } else if (strncmp(runargs[kk],"-dloginfo",8)==0){ 00861 iloginfo=atoi(runargs[kk+1]); 00862 } 00863 } 00864 info=DSDPLogInfoAllow(iloginfo,0); 00865 return 0; 00866 } 00867 00868 #undef __FUNCT__ 00869 #define __FUNCT__ "TCheckArgs" 00870 static int TCheckArgs(DSDP dsdp,SDPCone sdpcone, int m,int nargs,char *runargs[]){ 00871 00872 int i,kk, info; 00873 int printlevel=10; 00874 double *yy0; 00875 00876 for (kk=1; kk<nargs-1; kk++){ 00877 if (strncmp(runargs[kk],"-y0",7)==0){ 00878 DSDPCALLOC2(&yy0,double,m,&info); 00879 for (i=0;i<m;i++) yy0[i]=0.0; 00880 info = ReadInitialPoint(runargs[kk+1],m,yy0); 00881 for (i=0;i<m;i++){info = DSDPSetY0(dsdp,i+1,yy0[i]);} 00882 DSDPFREE(&yy0,&info); 00883 } else if (strncmp(runargs[kk],"-params",6)==0){ 00884 info=DSDPReadOptions(dsdp,runargs[kk+1]); 00885 } else if (strncmp(runargs[kk],"-print",6)==0){ 00886 printlevel=atoi(runargs[kk+1]); 00887 } 00888 } 00889 00890 info=DSDPSetOptions(dsdp,runargs,nargs); 00891 /* if (lpcone){info=LPConeScaleBarrier(lpcone,lpb); */ 00892 if (rank==0){info=DSDPSetStandardMonitor(dsdp,printlevel);} 00893 if (rank==0){info=DSDPSetFileMonitor(dsdp,printlevel);} 00894 return(0); 00895 } 00896 00897 #undef __FUNCT__ 00898 #define __FUNCT__ "ReadInitialPoint" 00899 static int ReadInitialPoint(char* filename, int m, double yy0[]) 00900 { 00901 FILE *fp; 00902 int i,count; 00903 00904 fp = fopen(filename,"r"); 00905 if(!fp) 00906 { printf("\n Cannot open file containing initial dual solution %s",filename); 00907 return 0; 00908 } 00909 00910 for(count=0,i=0;(i < m) &&(!feof(fp));i++){ 00911 if(fscanf(fp,"%lf",&yy0[i] ) ){ count++; } 00912 } 00913 00914 if (count < m){ 00915 printf("Warning: reading initial y vector: \n"); 00916 printf(" Expecting %d entries but read only %d entries \n",m,count); 00917 } 00918 fclose(fp); 00919 return 0; 00920 } /* ReadUserY */