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++){
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){
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
00158 for (i=0;i<m;i++){info = DSDPSetDualObjective(dsdp,i+1,dddd.dobj[i]);}
00159
00160
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
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 ){
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 ){
00192 info=SDPConeSetADenseVecMat(sdpcone,j,i,n,1.0,dddd.nnz+spot,ijnnz);CHKDATA(i,j,info);
00193 } else {
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
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
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
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
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 }
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
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
00426 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){
00427 fgets(thisline,BUFFERSIZ,fp); line++;
00428 }
00429
00430 if (sscanf(thisline,"%d",&m)<1){
00431 printf("Error: line %d. Number of constraints not given.\n",line);
00432 return(1);
00433 }
00434
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
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
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
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) ){
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
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
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
00845
00846
00847 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
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 }