DSDP
readsdpa.c
Go to the documentation of this file.
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 */