1 #include "dsdp5.h"
2 #include <string.h>
3 #include <stdio.h>
4 #include <math.h>
5 #include <stdlib.h>
6 #include "pdsdp5scalapack.h"
7
8 /*! \file readsdpa.c
9 \brief Read SDPA data files, pass data into DSDP solver, and print solution.
10 */
11
12 /*
13 static char help[]="\n";
14 */
15 static char help[]="\
16 DSDP Usage: dsdp5 filename <sdpa format file> \n\
17 -print <10> - print information at each k iteration\n\
18 -save <Solution File in SDPA format> \n\
19 -fout <filename> to print standard monitor to a file\n\
20 -y0 <initial solution file> \n\
21 -benchmark <file containing names of SDPA files> \n\
22 -directory <directory containing benchmark SDPA files> \n\
23 -suffix <add to each benchmark problem name> \n\
24 -dloginfo <0> - print more information for higher numbers \n\
25 -dlogsummary <1> - print timing information \n\
26 -help for this help message\n";
27
28 static int qusort(int[],int[],int[],double[],int,int);
29 static int partition(int[],int[],int[],double[],int, int);
30 static int Parseline(char *,int *,int *,int *,int *,double *, int *);
31 static int ReadInitialPoint(char*, int, double[]);
32 static int TCheckArgs0(DSDP,SDPCone,int,int,char *[]);
33 static int TCheckArgs(DSDP,SDPCone,int,int,char *[]);
34 static int CheckForConstantMat(double[],int, int);
35 static int CountNonzeroMatrices(int, int[],int[], int*);
36
37 typedef struct{
38 char sformat;
39 int blocksize;
40 } DBlock;
41
42 typedef struct{
43 int *block,*constraint,*matind;
44 double*nnz;
45 char *sformat;
46 int totalnonzeros;
47 double *dobj,*y0;
48 char *conetypes;
49 int *blocksizes;
50 int m; int n; int nblocks;
51 int lpn,lpspot,lpblock,lpnnz;
52 int *lpi,*lui,*cmap;
53 double cnorm;
54 double fixedvari;
55 double fixedvard,xout;
56 } DSDPData;
57
58 static int ReadSDPA2(char*,DSDPData*);
59 static int GetMarkers(int, int, int*, int*, int*);
60 static int ComputeY0(DSDP,DSDPData);
61 static int rank=0;
62 int ReadSDPAFile(int argc,char *argv[]);
63
64 #define CHKDATA(a,b,c) { if (c){ printf("Possible problem in variable %d, block %d. \n",a+1,b+1);} }
65
66
67 #undef __FUNCT__
68 #define __FUNCT__ "main"
main(int argc,char * argv[])69 int main(int argc,char *argv[]){
70 int info;
71 MPI_Init(&argc,&argv);
72 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
73 DSDPSetRank(rank);
74 info=ReadSDPAFile(argc,argv);
75 MPI_Finalize( );
76 return info;
77 }
78
79 #undef __FUNCT__
80 #define __FUNCT__ "ReadSDPAFile"
81 /*!
82 \fn int ReadSDPAFile(int argc,char *argv[]);
83 \brief Read SDPA formatted file and solve the semidefinite program.
84 \param argc number of command line arguments
85 \param argv command line arguments
86 \ingroup Examples
87 */
ReadSDPAFile(int argc,char * argv[])88 int ReadSDPAFile(int argc,char *argv[]){
89
90 int i,j,m,n,np,its,info;
91 int spot,ijnnz,nzmats,sdpnmax,sdpn,stat1,printsummary=1;
92 int runbenchmark=0,saveit=0,justone=1,fileout=0;
93 double t1,t2,t3,t4,t5,dd,yhigh;
94 double derr[6],dnorm[3];
95 double ddobj,ppobj;
96 char problemname[100],thisline[100], filename[300],savefile[100];
97 char directory[100]="/home/benson/sdpexamples/sdplib/";
98 char outputfile[50]="",suffix[20]=".dat-s", tablename[20]="results-dsdp-5.7";
99 char success='f',sformat;
100 FILE *fp1=0,*fp2=0,*fout;
101 DSDPData dddd;
102 DSDP dsdp;
103 DSDPTerminationReason reason;
104 DSDPSolutionType pdfeasible;
105 SDPCone sdpcone=0;
106 LPCone lpcone=0;
107 int *ittt,sspot;
108
109 if (argc<2){ printf("%s",help); DSDPPrintOptions(); return(0); }
110 for (i=0; i<argc; i++){ if (strncmp(argv[i],"-help",5)==0){
111 printf("%s",help); DSDPPrintOptions(); return(0);}}
112 for (i=1; i<argc; i++){ printf("%s ",argv[i]); } printf("\n");
113 for (i=1; i<argc-1; i++){ /* Are we reading a file with a lot of problem names? */
114 if (strncmp(argv[i],"-benchmark",8)==0){
115 strncpy(thisline,argv[i+1],90); fp1=fopen(thisline,"r");runbenchmark=1; justone=0;
116 };
117 if (strncmp(argv[i],"-directory",8)==0){strncpy(directory,argv[i+1],90);}
118 if (strncmp(argv[i],"-table",4)==0){strncpy(tablename,argv[i+1],90);};
119 if (strncmp(argv[i],"-suffix",4)==0){strncpy(suffix,argv[i+1],20);};
120 if (strncmp(argv[i],"-save",5)==0){ strncpy(savefile,argv[i+1],40);saveit=1;};
121 if (strncmp(argv[i],"-dlogsummary",8)==0){printsummary=atoi(argv[i+1]);}
122 if (rank==0&&strncmp(argv[i],"-fout",5)==0){ strncpy(outputfile,argv[i+1],45);fileout=1;};
123 }
124
125
126 while ((runbenchmark && !feof(fp1)) || justone==1){
127 justone=0;
128 if (runbenchmark){ /* Get filename with the data */
129 fgets(thisline,100,fp1); if (sscanf(thisline,"%s",problemname)<1) break;
130 strncpy(filename,directory,90); strncat(filename,problemname,90);strncat(filename,suffix,90);
131 printf("%s\n",problemname);
132 } else {
133 strncpy(filename,argv[1],90);
134 strncpy(problemname,argv[1],90);
135 }
136
137 if (rank==0 && fileout){
138 dsdpoutputfile=fopen(outputfile,"a");
139 fprintf(dsdpoutputfile,"%s\n",problemname);
140 } else {dsdpoutputfile=0;fileout=0;}
141 DSDPTime(&t1);
142
143 info=ReadSDPA2(filename, &dddd); m=dddd.m;
144 if (info){ printf("Problem reading SDPA file\n"); return 1;}
145 DSDPTime(&t2);
146 if (rank==0){
147 printf("\nVariables %d \n",dddd.m);
148 printf("Matrix Blocks: %d, ",dddd.nblocks);
149 printf("Total Number of Constraints: %d \n",dddd.n);
150 printf("Nonzeros in Constraints: %d\n\n",dddd.totalnonzeros);
151 printf("Read Data File into Buffer: %4.3e seconds\n",t2-t1);
152 }
153
154 for (i=0; i<argc-1; i++){
155 if (strncmp(argv[i],"-dloginfo",8)==0){ info=DSDPLogInfoAllow(atoi(argv[i+1]),0);};
156 }
157
158 info = DSDPCreate(dddd.m,&dsdp);
159 info = DSDPCreateSDPCone(dsdp,dddd.nblocks,&sdpcone);
160 /* Set Dual objective vector */
161 for (i=0;i<m;i++){info = DSDPSetDualObjective(dsdp,i+1,dddd.dobj[i]);}
162
163 /* Set initial point */
164 for (i=0; i<m; i++)
165 if (dddd.dobj[i]> 0.0) dddd.y0[i]=-0.0; else dddd.y0[i]=0.0;
166 for (i=0; i<m; i++){info = DSDPSetY0(dsdp,i+1,dddd.y0[i]);}
167 info=ComputeY0(dsdp,dddd);
168 if (dddd.fixedvari){
169 info = DSDPSetY0(dsdp,(int)dddd.fixedvari,dddd.fixedvard);
170 printf("Fixed: %2.0f %4.2f ?\n",dddd.fixedvari,dddd.fixedvard);
171 DSDPSetFixedVariables(dsdp,&dddd.fixedvari,&dddd.fixedvard,&dddd.xout,1);
172 }
173
174 spot=0;ijnnz=0;np=0;sdpnmax=1;sdpn=0;stat1=1;
175 /* Insert the SDP data */
176 for (j=0;j<dddd.nblocks; j++){
177 if (dddd.conetypes[j]=='S'){
178 n=dddd.blocksizes[j];
179 sformat=dddd.sformat[j];
180 info=CountNonzeroMatrices(j+1,dddd.block+spot,dddd.constraint+spot,&nzmats);
181 info=SDPConeSetBlockSize(sdpcone,j,n); DSDPCHKERR(info);
182 info=SDPConeSetSparsity(sdpcone,j,nzmats); DSDPCHKERR(info);
183 info=SDPConeSetStorageFormat(sdpcone,j,sformat); DSDPCHKERR(info);
184 np+=n; sdpn+=n;
185 if (sdpnmax<n) sdpnmax=n;
186 if (stat1<nzmats) stat1=nzmats;
187 for (i=0; i<=m; i++){
188 info=GetMarkers(j+1,i,dddd.block+spot,dddd.constraint+spot,&ijnnz);
189 if (0==1){
190 } else if ( ijnnz==0 ){ /* info=DSDPSetZeroMat(dsdp,j,i,n); */
191 } else if (CheckForConstantMat(dddd.nnz+spot,ijnnz,n)){
192 info=SDPConeSetConstantMat(sdpcone,j,i,n,dddd.nnz[spot+1]);CHKDATA(i,j,info);
193 if(sformat=='P'){info=SDPConeSetXArray(sdpcone,j,n,dddd.nnz+spot,n*(n+1)/2);}
194 } else if (sformat=='P' && ijnnz==n*(n+1)/2 ){ /* check for dense matrix */
195 info=SDPConeSetADenseVecMat(sdpcone,j,i,n,1.0,dddd.nnz+spot,ijnnz);CHKDATA(i,j,info);
196 } else { /* sparse matrix */
197 info=SDPConeSetASparseVecMat(sdpcone,j,i,n,1.0,0,dddd.matind+spot,dddd.nnz+spot,ijnnz);CHKDATA(i,j,info);
198 }
199 if (0==1){ info=SDPConeViewDataMatrix(sdpcone,j,i);}
200 spot+=ijnnz;
201 /* SDPConeScaleBarrier(sdpcone,j,j+1.0); */
202 }
203 if (0==1){info=SDPConeView2(sdpcone);}
204 } else if (dddd.conetypes[j]=='L'){
205 info=DSDPCreateLPCone(dsdp,&lpcone); sformat='P';
206 info=SDPConeSetStorageFormat(sdpcone,j,sformat); DSDPCHKERR(info);
207 n=dddd.blocksizes[j];
208 np+=n;
209 sspot=spot;
210 DSDPCALLOC2(&ittt,int,(m+2),&info);
211 for (i=0;i<=m;i++){ittt[i]=0;}
212 for (i=0;i<=m;i++){
213 info=GetMarkers(j+1,i,dddd.block+spot,dddd.constraint+spot,&ijnnz);
214 ittt[i+1]=ijnnz; spot+=ijnnz;
215 }
216 for (i=1;i<=m;i++)ittt[i+1]+=ittt[i];
217 info=LPConeSetData(lpcone,n,ittt,dddd.matind+sspot,dddd.nnz+sspot);CHKDATA(i,0,info);
218 if (0==1){info=LPConeView(lpcone);}
219 if (0==1){info=LPConeView2(lpcone);}
220 /* info=DSDPFree(&ittt); */
221 }
222 }
223 if (0==1){
224 BCone bcone;
225 info=DSDPCreateBCone(dsdp, &bcone);
226 info=BConeAllocateBounds(bcone,2*m);
227 for (i=0;i<m;i++){
228 info=BConeSetUpperBound(bcone,i+1,10);
229 }
230 for (i=0;i<m;i++){
231 info=BConeSetLowerBound(bcone,i+1,-10);
232 }
233 }
234
235 DSDPTime(&t3);
236 if (rank==0){printf("DSDP Set Data: %4.3e seconds\n",t3-t2);}
237
238 its=(m-2)/sdpnmax;
239 if (np<100 && its==0) its=1;
240 if (dddd.lpn>=m && its==0) its=1;
241 if (its>=1) its++;
242 its=its*its;
243 if (m<2000 && its>10) its=10;
244 if (its>12) its=12;
245
246 info=DSDPReuseMatrix(dsdp,its);
247
248 DSDPFREE(&dddd.blocksizes,&info);
249 DSDPFREE(&dddd.sformat,&info);
250 DSDPFREE(&dddd.dobj,&info);
251 DSDPFREE(&dddd.y0,&info);
252 DSDPFREE(&dddd.conetypes,&info);
253 DSDPFREE(&dddd.constraint,&info);
254 DSDPFREE(&dddd.block,&info);
255
256 info=DSDPGetDataNorms(dsdp, dnorm);
257 if (dnorm[0]==0){
258 info=DSDPSetR0(dsdp,np);
259 info=DSDPSetGapTolerance(dsdp,1e-3);
260 info=DSDPSetYBounds(dsdp,-1.0,1.0);
261 } else {
262 }
263 info = TCheckArgs0(dsdp,sdpcone,dddd.m,argc,argv);
264 info = TCheckArgs(dsdp,sdpcone,dddd.m,argc,argv);
265
266 info=PDSDPUseSCALAPACKLinearSolver(dsdp);
267 info = DSDPSetup(dsdp); if (info){ printf("\nProblem Setting problem. Likely insufficient memory\n"); return 1;}
268 if (0==1){info=SDPConeCheckData(sdpcone);}
269
270
271 DSDPTime(&t4);
272 if (rank==0){
273 printf("DSDP Process Data: %4.3e seconds\n\n",t4-t3);
274 printf("Data Norms: C: %4.2e, A: %4.2e, b: %4.2e\n",dnorm[0],dnorm[1],dnorm[2]);
275 info=DSDPGetScale(dsdp,&ddobj);
276 printf("Scale C: %4.2e\n\n",ddobj);
277 info=DSDPGetPotentialParameter(dsdp,&ddobj);
278 printf("Potential Parameter: %4.2f\n",ddobj);
279 info=DSDPGetReuseMatrix(dsdp,&its);
280 printf("Reapply Schur matrix: %d\n\n",its);
281 }
282 if (0==1){info=DSDPPrintData(dsdp,sdpcone,lpcone);}
283
284 info = DSDPSolve(dsdp);
285 if (info){ printf("\nNumerical errors encountered in DSDPSolve(). \n");}
286
287 info=DSDPStopReason(dsdp,&reason);
288 if (reason!=DSDP_INFEASIBLE_START){
289 info=DSDPComputeX(dsdp);DSDPCHKERR(info);
290 }
291 info=DSDPStopReason(dsdp,&reason);
292 info=DSDPGetSolutionType(dsdp,&pdfeasible);
293
294 DSDPTime(&t5);
295
296 if (rank==0){
297
298 if (reason == DSDP_CONVERGED){
299 printf("DSDP Converged. \n");
300 success='s';
301 } else if ( reason == DSDP_UPPERBOUND ){
302 printf("DSDP Terminated Because Dual Objective Exceeded its Bound\n");
303 success='s';
304 } else if ( reason == DSDP_SMALL_STEPS ){
305 printf("DSDP Terminated Due to Small Steps\n");
306 success='f';
307 } else if ( reason == DSDP_MAX_IT){
308 printf("DSDP Terminated Due Maximum Number of Iterations\n");
309 success='f';
310 } else if ( reason == DSDP_INFEASIBLE_START){
311 printf("DSDP Terminated Due to Infeasible Starting Point\n");
312 success='f';
313 } else if ( reason == DSDP_INDEFINITE_SCHUR_MATRIX){
314 printf("DSDP Terminated Due to Indefinite Schur Complement\n");
315 success='f';
316 } else {
317 printf("DSDP Finished\n");
318 success='f';
319 }
320
321 if (pdfeasible == DSDP_UNBOUNDED ){
322 printf("DSDP Dual Unbounded, Primal Infeasible\n");
323 } else if ( pdfeasible == DSDP_INFEASIBLE ){
324 printf("DSDP Primal Unbounded, Dual Infeasible\n");
325 }
326
327
328 info=DSDPGetDObjective(dsdp,&ddobj);
329 info=DSDPGetPObjective(dsdp,&ppobj);
330 printf("\nP Objective : %16.8e \n",ppobj);
331 printf("DSDP Solution: %16.8e \n\n",ddobj);
332 printf("DSDP Solve Time: %4.3e seconds\n",t5-t4);
333 printf("DSDP Preparation and Solve Time: %4.3e seconds\n\n",t5-t3);
334
335 info=DSDPGetFinalErrors(dsdp,derr);
336 info=DSDPGetIts(dsdp,&its);
337 if (1 || runbenchmark){
338 fp2=fopen(tablename,"a");
339 if (pdfeasible==DSDP_UNBOUNDED){
340 fprintf(fp2," %-18s & %4d & %4d & infeasible & unbounded & %4.0e & %c & %3d & %6.2f \\\\\n",problemname,m,np,derr[0],success,its,t5-t3);
341 } else if (pdfeasible==DSDP_INFEASIBLE){
342 fprintf(fp2," %-18s & %4d & %4d & unbounded & infeasible & %4.0e & %c & %3d & %6.2f \\\\\n",problemname,m,np,derr[0],success,its,t5-t3);
343 } else {
344 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);
345 }
346 fclose(fp2);
347 }
348
349 /* info=DSDPComputeMinimumXEigenvalue(dsdp,&derr[1]); */
350 printf("\nP Infeasible: %8.2e \n",derr[0]);
351 printf("D Infeasible: %8.2e \n",derr[2]);
352 printf("Minimal P Eigenvalue: %6.2e \n",derr[1]);
353 printf("Minimal D Eigenvalue: 0.00 \n");
354 printf("Relative P - D Objective values: %4.2e \n",derr[4]);
355 printf("Relative X Dot S: %4.2e \n",derr[5]);
356 info=DSDPGetYBounds(dsdp,&dd,&yhigh);
357 info=DSDPGetYMaxNorm(dsdp,&dd);
358 printf("\nMax Y: %10.8e, Bounded by %6.1e\n",dd,yhigh);
359 info=DSDPGetTraceX(dsdp,&dd);
360 printf("Trace X: %4.8e, ",dd);
361 info=DSDPGetPenaltyParameter(dsdp,&dd);
362 printf("Bounded by Penalty Parameter: %4.1e \n\n",dd);
363
364 if (printsummary){ DSDPEventLogSummary();}
365 printf("--- DSDP Finished ---\n\n");
366
367 if (saveit){
368 fout=fopen(savefile,"w");
369 /* fprintf(fout,"** %s \n",filename); Deleted */
370 info= DSDPPrintSolution(fout,dsdp,sdpcone,lpcone);
371 if (dddd.fixedvari){
372 sspot=dddd.nblocks+1,dd=dddd.xout;
373 fprintf(fout,"1 %d 1 1 1.0e-11\n1 %d 2 2 1.0e-11\n",sspot,sspot);
374 fprintf(fout,"2 %d 1 1 %12.8e\n",sspot,DSDPMax(1.0e-10,dd));
375 fprintf(fout,"2 %d 2 2 %12.8e\n",sspot,DSDPMax(1e-10,-dd));
376 }
377 fclose(fout);
378 }
379 }
380
381 info = DSDPDestroy(dsdp);
382 DSDPFREE(&dddd.matind,&info);
383 DSDPFREE(&dddd.nnz,&info);
384
385 if (fileout){fclose(dsdpoutputfile);}
386 if (0){ DSDPMemoryLog();}
387
388 }
389 if (runbenchmark){ fclose(fp1);}
390
391 return 0;
392 } /* main */
393
394
395
396 #define BUFFERSIZ 4000
397 #undef __FUNCT__
398 #define __FUNCT__ "ReadSDPA2"
ReadSDPA2(char * filename,DSDPData * ddd)399 static int ReadSDPA2(char *filename, DSDPData*ddd){
400 FILE*fp;
401 char ctmp,refline[BUFFERSIZ]="*",thisline[BUFFERSIZ]="*";
402 int info,tline,line=0;
403 int i,k,m,n;
404 /* int spot,nzmark, */
405 int bigint=1000000;
406 int i1,nblk,nmat,col,row;
407 int np=0,nblocks;
408 int nargs,nonzero;
409 double val;
410
411 fp=fopen(filename,"r");
412 if (!fp){
413 printf("Cannot open file %s !",filename); return(1);
414 }
415
416 /* Read comments */
417 while(!feof(fp) && (thisline[0] == '*' || thisline[0] == '"') ){
418 fgets(thisline,BUFFERSIZ,fp); line++;
419 }
420 /* Read number of constraints */
421 if (sscanf(thisline,"%d",&m)<1){
422 printf("Error: line %d. Number of constraints not given.\n",line);
423 return(1);
424 }
425 /* Read number of blocks */
426 fgets(thisline,BUFFERSIZ,fp); line++;
427 if (sscanf(thisline,"%d",&nblocks)!=1){
428 printf("Error: line %d. Number of blocks not given.\n",line);
429 return(1);
430 }
431 ddd->lpn=0;ddd->lpspot=0;ddd->lpblock=0;ddd->cnorm=0;
432 /* Read block sizes */
433 DSDPCALLOC2(&ddd->sformat,char, (nblocks+1),&info);
434 DSDPCALLOC2(&ddd->blocksizes,int, (nblocks+1),&info);
435 DSDPCALLOC2(&ddd->conetypes,char, (nblocks+1),&info );
436 line++;
437 for (i=0;i<nblocks; i++){
438 if (fscanf(fp,"{")==1 || fscanf(fp,"(")==1 || fscanf(fp,",")==1 ){
439 i--;
440 } else if (fscanf(fp,"%d",&col)==1){
441 if (col>0) { ddd->blocksizes[i]=col; np+=col; ddd->conetypes[i]='S';
442 } else if (col>0){ddd->blocksizes[i]=-col; np+=-col; ddd->conetypes[i]='S';
443 } else if (col<0){ddd->blocksizes[i]=-col; np += -col; ddd->conetypes[i]='L';ddd->lpn=-col;ddd->lpblock=i;
444 } else { ddd->blocksizes[i]=0; ddd->conetypes[i]='N';}
445 if (ddd->blocksizes[i]<10){ddd->sformat[i]='U';} else {ddd->sformat[i]='P';}
446 }
447 else{ printf("Error block sizes, line %d",line); return(1);}
448 }
449 if (ddd->blocksizes[nblocks-1]==0) nblocks--;
450 fgets(thisline,BUFFERSIZ,fp);
451
452 /* Read objective vector */
453 DSDPCALLOC2(&ddd->y0,double,m,&info);
454 DSDPCALLOC2(&ddd->dobj,double,m,&info);
455 line++;
456 for (i=0;i<m;i++){
457 if (fscanf(fp,",")==1){
458 i--;
459 continue;
460 }
461 while (fscanf(fp,"%lg",&val)!=1){
462 fscanf(fp,"%c",&ctmp);
463 if (ctmp=='\n'){
464 printf("Constraints: %d, Blocks: %d\n",m,nblocks);
465 printf("Error reading objective, line %d, i=%d \n",line,i); return 1;
466 }
467 }
468 ddd->dobj[i]=val;
469 }
470 fgets(thisline,BUFFERSIZ,fp);
471 tline=line;
472
473 nargs=5; nonzero=0;
474 while(!feof(fp)){
475 thisline[0]='\0';
476 nmat=-1; nblk=-1; row=-1; col=-1; val=0.0;
477 fgets(thisline,BUFFERSIZ,fp); line++;
478 info = Parseline(thisline,&nmat,&nblk,&row,&col,&val,&nargs);
479 if (!feof(fp)&&nargs!=5&&nargs>0){
480 printf("Error: line %d \n%s\n",line,thisline);return 1;}
481 if (nargs==5 && val!=0.0){
482 nonzero++;
483 i1=row*(row+1)/2 + col;
484 if (row >= ddd->blocksizes[nblk-1] || col >= ddd->blocksizes[nblk-1] ) {
485 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);
486 return 1;
487 }
488 if (row<0 || col<0){
489 printf("Data Error in line: %d. Row %d or col %d <= 0 \n%s",line,row+1,col+1,thisline);
490 return 1;
491 }
492 if (nmat>m || nmat<0){
493 printf("Data Error in line: %d. Is Var 0 <= %d <= %d \n%s",line,nmat,m,thisline);
494 return 1;
495 }
496 if (nblk>nblocks || nblk<0){
497 printf("Data Error in line: %d. Is Block 0 <= %d <= %d \n%s",line,nmat,m,thisline);
498 return 1;
499 }
500 } else if (nargs==5 && val==0.0){
501 } else if (nargs==0){
502 } else {
503 printf("Problem Reading SDPA file at line %d: %s\n",line, thisline);
504 }
505 }
506
507 /* Allocate memory for the data */
508 nonzero++;
509 DSDPCALLOC2(&ddd->matind,int,nonzero,&info);
510 DSDPCALLOC2(&ddd->nnz,double,nonzero,&info);
511 DSDPCALLOC2(&ddd->block,int,nonzero,&info);
512 DSDPCALLOC2(&ddd->constraint,int,nonzero,&info);
513 nonzero--;
514
515 fseek(fp,0,SEEK_SET);
516 line=0;
517 for (i=0;i<tline;i++){ctmp='*'; while (ctmp!='\n') fscanf(fp,"%c",&ctmp); line++;}
518
519 nargs=5;k=0;
520 while(!feof(fp) /* && nargs==5 */){
521 thisline[0]='\0';
522 fgets(thisline,BUFFERSIZ,fp);
523 if (k==0){strncpy(refline,thisline,BUFFERSIZ-1); }
524 info = Parseline(thisline,&nmat,&nblk,&row,&col,&val,&nargs);
525 if (!feof(fp)&&nargs!=5&&nargs<0){
526 /* if (k>=nonzero && !feof(fp) ){ */
527 printf("Problem Reading SDPA file at line %d: %s\n",line, thisline);
528 printf("Problem could be earlier in file \n");
529 printf("The first recorded matix nonzero in the file occured at line %d: \n %s",tline,refline);
530 printf(" Please check data file\n");
531 return 1;
532 }
533 if (nargs==5 && val!=0.0){
534 if (row>col){
535 printf("Warning: Line: %d Row < Column. %s \n",line,thisline);
536 }
537 i=row;row=col;col=i;
538 n=ddd->blocksizes[nblk-1];
539 if (nmat==0) {val=-val;}
540 if (ddd->conetypes[nblk-1]=='S'){
541 /* if (row==col) val/=2; */
542 ddd->matind[k]=row*(row+1)/2 + col;
543 if (ddd->sformat[nblk-1]=='U'){ddd->matind[k]=row*n + col;}
544 } else {
545 ddd->matind[k]=col;
546 }
547 ddd->block[k]=nblk;
548 ddd->constraint[k]=nmat;
549 ddd->nnz[k]=val;
550 k++;
551 } else if (nargs==5 && val==0.0){
552 } else if (nargs==0){
553 } else {
554 printf("Problem Reading SDPA file at line %d: %s\n",line, thisline);
555 printf("Problem could be earlier in file \n");
556 printf("The first recorded matix nonzero in the file occured at line %d: \n %s",tline,refline);
557 printf(" Please check data file\n");
558 return 1;
559 }
560 }
561 ddd->block[k]=nblocks+1; ddd->constraint[k]=m+2;
562 ddd->matind[k]=10000000; ddd->nnz[k]=0.0;
563
564 qusort(ddd->block,ddd->constraint,ddd->matind,ddd->nnz,0,nonzero-1);
565
566 for (i=0;i<nonzero-1; i++){
567 while (i<nonzero-1 && ddd->matind[i]==ddd->matind[i+1] &&
568 ddd->constraint[i]==ddd->constraint[i+1] &&
569 ddd->block[i]==ddd->block[i+1] ){
570 printf("DSDPError: Reading Input File:\n");
571 printf("Possible problem with data input file: Double Entry: \n");
572 printf(" %d %d %d %2.10e\n",
573 ddd->constraint[i],ddd->block[i],ddd->matind[i]+1,ddd->nnz[i]);
574 printf(" %d %d %d %2.10e\n\n",
575 ddd->constraint[i+1],ddd->block[i+1],ddd->matind[i+1]+1,ddd->nnz[i+1]);
576 for (k=i+1;k<nonzero-1;k++){
577 ddd->constraint[k]=ddd->constraint[k+1]; ddd->block[k]=ddd->block[k+1];
578 ddd->matind[k]=ddd->matind[k+1];ddd->nnz[k]=ddd->nnz[k+1];
579 }
580 ddd->constraint[nonzero-1]=bigint;ddd->nnz[nonzero-1]=0;
581 nonzero--;
582 }
583 }
584
585 ddd->fixedvari=0;ddd->fixedvard=0;
586 if (ddd->lpblock>0){
587 int spot;
588 if (ddd->blocksizes[ddd->lpblock]==2){
589 i=0;k=0;
590 while (ddd->block[i]<=ddd->lpblock && i<nonzero){ i++;} spot=i;
591 while (ddd->block[i]==ddd->lpblock+1 && i<nonzero){ i++;k++;}
592 if (k==4){
593 if (ddd->constraint[spot]==ddd->constraint[spot+1] &&
594 ddd->constraint[spot+2]==ddd->constraint[spot+3] &&
595 ddd->matind[spot]==ddd->matind[spot+2] &&
596 ddd->matind[spot+1]==ddd->matind[spot+3] &&
597 fabs(ddd->nnz[spot+2])==1.0 && fabs(ddd->nnz[spot+3])==1.0 &&
598 fabs(ddd->nnz[spot] + ddd->nnz[spot+1]) <=1e-6 ){
599 ddd->fixedvari=ddd->constraint[spot+2];
600 ddd->fixedvard=ddd->nnz[spot]/ddd->nnz[spot+2];
601 nblocks--;ddd->lpblock=0;
602 }
603 }
604 }
605 }
606
607 ddd->totalnonzeros=nonzero;
608 for (ddd->n=0,i=0;i<nblocks;i++) ddd->n += ddd->blocksizes[i];
609 ddd->m=m;
610 ddd->nblocks=nblocks;
611 fclose(fp);
612 return 0;
613 }
614
615
616 #undef __FUNCT__
617 #define __FUNCT__ "Parseline"
Parseline(char thisline[],int * nmat,int * nblk,int * row,int * col,double * value,int * nargs)618 static int Parseline(char thisline[],int *nmat,int *nblk,int *row,
619 int *col,double *value, int *nargs){
620
621 int temp;
622 int rtmp,coltmp;
623
624 *nmat=-1;*nblk=-1;rtmp=-1;coltmp=-1;*value=0.0;
625 temp=sscanf(thisline,"%d %d %d %d %lg",nmat,nblk,&rtmp,&coltmp,value);
626 if (temp==5) *nargs=5;
627 else *nargs=0;
628 *row=rtmp-1; *col=coltmp-1;
629
630 return(0);
631 }
632
633
partition(int list1[],int list2[],int list3[],double list5[],int lstart,int lend)634 static int partition(int list1[], int list2[], int list3[], double list5[], int lstart, int lend){
635 int k=lend;
636 int pivot1=list1[k],pivot2=list2[k],pivot3=list3[k];
637 double pivot5 = list5[k];
638 int bottom = lstart-1, top = lend;
639 int done = 0;
640 int ordered=1;
641 while (!done){
642
643 while (!done) {
644 bottom = bottom+1;
645
646 if (bottom == top){
647 done = 1;
648 break;
649 }
650 if ( list1[bottom] > pivot1 ||
651 (list1[bottom] == pivot1 && list2[bottom] > pivot2) ||
652 (list1[bottom] == pivot1 && list2[bottom] == pivot2 &&
653 list3[bottom] > pivot3) ){
654 list1[top] = list1[bottom];
655 list2[top] = list2[bottom];
656 list3[top] = list3[bottom];
657 list5[top] = list5[bottom];
658 ordered=0;
659 break;
660 }
661 }
662 while (!done){
663 top = top-1;
664
665 if (top == bottom){
666 done = 1;
667 break;
668 }
669 if ( list1[top] < pivot1 ||
670 (list1[top] == pivot1 && list2[top] < pivot2) ||
671 (list1[top] == pivot1 && list2[top] == pivot2 && list3[top] < pivot3)){
672 list1[bottom] = list1[top];
673 list2[bottom] = list2[top];
674 list3[bottom] = list3[top];
675 list5[bottom] = list5[top];
676 ordered=0;
677 break;
678 }
679 }
680 }
681 list1[top] = pivot1;
682 list2[top] = pivot2;
683 list3[top] = pivot3;
684 list5[top] = pivot5;
685
686 ordered=0;
687
688 if (bottom==lend){
689 ordered=1;
690 for (k=lstart;k<lend-1;k++){
691 if ( (list1[k] > list1[k+1]) ||
692 (list1[k] == list1[k+1] && list2[k] > list2[k+1]) ||
693 (list1[k] == list1[k+1] && list2[k] == list2[k+1] && list3[k] > list3[k+1]) ){
694 ordered=0;
695 break;
696 }
697 }
698 }
699 if (ordered && lend-lstart>2){
700 top=(lend+lstart)/2;
701 }
702 return top;
703 }
704
705
706
qusort(int list1[],int list2[],int list3[],double list5[],int lstart,int lend)707 static int qusort(int list1[], int list2[], int list3[], double list5[], int lstart, int lend){
708 int split;
709 if (lstart < lend){
710 split = partition(list1, list2, list3, list5,lstart, lend);
711 qusort(list1, list2, list3, list5, lstart, split-1);
712 qusort(list1, list2, list3, list5, split+1, lend);
713 } else {
714 return 0;
715 }
716 return 0;
717 }
718
719
720
721 #undef __FUNCT__
722 #define __FUNCT__ "GetMarkers"
GetMarkers(int block,int constraint,int blockn[],int constraintn[],int * m3)723 static int GetMarkers(int block, int constraint, int blockn[], int constraintn[],
724 int*m3){
725 int i=0;
726 while (blockn[i]==block && constraintn[i]==constraint){ i++;}
727 *m3=i;
728 return 0;
729 }
730
731 #undef __FUNCT__
732 #define __FUNCT__ "CountNonzeroMatrices"
CountNonzeroMatrices(int block,int blockn[],int constraintn[],int * m3)733 static int CountNonzeroMatrices(int block, int blockn[], int constraintn[], int*m3){
734 int i=0,cvar=-1,nnzmats=0;
735 while (blockn[i]==block){
736 if (constraintn[i]>cvar){ cvar=constraintn[i];nnzmats++;}
737 i++;
738 }
739 *m3=nnzmats;
740 return 0;
741 }
742
743 #undef __FUNCT__
744 #define __FUNCT__ "CheckForConstantMat"
CheckForConstantMat(double v[],int nnz,int n)745 static int CheckForConstantMat(double v[],int nnz, int n){
746 int i; double vv;
747 if (n<=1){ return 0; }
748 if (nnz!=(n*n+n)/2){ return 0; }
749 for (vv=v[0],i=1;i<nnz;i++){
750 if (v[i]!=vv){ return 0;}
751 }
752 return 1;
753 }
754
ComputeY0(DSDP dsdp,DSDPData dddd)755 static int ComputeY0(DSDP dsdp,DSDPData dddd){
756 int i,ii,info,ijnnz=0,spot=0,ddiag=0,diag=0,n=dddd.blocksizes[0],m=dddd.m;
757 double aa,bb=0,ddmax=0,dd=0,cnorm=0;
758 char sformat=dddd.sformat[0];
759 if (dddd.nblocks>1) return 0;
760 if (dddd.fixedvari) return 0;
761
762 info=GetMarkers(1,0,dddd.block+spot,dddd.constraint+spot,&ijnnz);
763 for (i=0;i<ijnnz;i++){if (cnorm<fabs(dddd.nnz[i])) cnorm=fabs(dddd.nnz[i]);}
764 spot+=ijnnz;
765
766 for (i=1;i<=m;i++,spot+=ijnnz){
767 info=GetMarkers(1,i,dddd.block+spot,dddd.constraint+spot,&ijnnz);
768 ddiag=0;
769 if (ijnnz==n){
770 dd=dddd.nnz[spot]; ddiag=1;bb=dddd.dobj[i-1];
771 for (ii=0;ii<n;ii++){
772 if (sformat=='U'){
773 if (dddd.matind[spot+ii] != ii*n+ii || dddd.nnz[spot+i]!=dd){ ddiag=0;}
774 } else if (sformat=='P'){
775 if (dddd.matind[spot+ii] != ii*(ii+1)/2+ii || dddd.nnz[spot+i]!=dd){ ddiag=0;}
776 }
777 }
778 }
779 if (ddiag){
780 info=DSDPSetY0(dsdp,i,-100*n*cnorm*dddd.dobj[i-1]);
781 info=DSDPSetR0(dsdp,0);
782 info=DSDPSetZBar(dsdp,100*dd/bb*cnorm);
783 info=DSDPSetPotentialParameter(dsdp,5.0);
784 }
785 }
786
787 if (m==n){
788 spot=0; info=GetMarkers(1,0,dddd.block+spot,dddd.constraint+spot,&ijnnz); spot+=ijnnz;
789 dd=dddd.nnz[spot]; bb=dddd.dobj[0];
790 for (diag=1,i=0;i<n;i++,spot+=ijnnz){
791 info=GetMarkers(1,i+1,dddd.block+spot,dddd.constraint+spot,&ijnnz);
792 dd=dddd.nnz[spot];bb=dddd.dobj[i];
793 if (ddmax<bb/dd) ddmax=bb/dd;
794 if (sformat=='P'){
795 if (ijnnz!=1 || bb!=dddd.dobj[i] || dd!=dddd.nnz[spot] || dddd.matind[spot]!= ((i+1)*(i+2))/2-1) { diag=0; }
796 } else if (sformat=='U'){
797 if (ijnnz!=1 || bb!=dddd.dobj[i] || dd!=dddd.nnz[spot] || dddd.matind[spot]!= ((i)*n)+i) { diag=0; }
798 }
799 }
800 if (diag && cnorm*sqrt(1.0*n)<1e6){
801 for (i=0;i<n;i++){info = DSDPSetY0(dsdp,i+1,-10*sqrt(1.0*n)*cnorm);}
802 info=DSDPSetR0(dsdp,0); info=DSDPSetZBar(dsdp,100*ddmax*n*cnorm); info=DSDPSetPotentialParameter(dsdp,5.0); info=DSDPReuseMatrix(dsdp,0);
803 }
804 }
805 if (m==n+1){
806 spot=0; info=GetMarkers(1,0,dddd.block+spot,dddd.constraint+spot,&ijnnz); spot+=ijnnz;
807 dd=dddd.nnz[spot]; bb=dddd.dobj[0];diag=1;
808 info=GetMarkers(1,1,dddd.block+spot,dddd.constraint+spot,&ijnnz);
809 if (CheckForConstantMat(dddd.nnz+spot,ijnnz,n)){aa=dddd.nnz[0]; spot+=ijnnz;ii=2;} else {ii=1;}
810 for (i=0;i<n;i++,spot+=ijnnz){
811 info=GetMarkers(1,i+ii,dddd.block+spot,dddd.constraint+spot,&ijnnz);
812 dd=dddd.nnz[spot];bb=dddd.dobj[i+ii-1];
813 if (ddmax<bb/dd) ddmax=bb/dd;
814 if (sformat=='U'){
815 if (ijnnz!=1 || bb!=dddd.dobj[ii-1] || dd!=dddd.nnz[spot] || dddd.matind[spot]!= ((i)*(n))+i ) { diag=0; }
816 } else if (sformat=='P'){
817 if (ijnnz!=1 || bb!=dddd.dobj[ii-1] || dd!=dddd.nnz[spot] || dddd.matind[spot]!= ((i+1)*(i+2))/2-1) { diag=0; }
818 }
819 }
820 if (ii==1 && diag==1){
821 info=GetMarkers(1,m,dddd.block+spot,dddd.constraint+spot,&ijnnz);
822 if (CheckForConstantMat(dddd.nnz+spot,ijnnz,n)){aa=dddd.nnz[spot];} else {diag=0;}
823 }
824 if (diag && cnorm*sqrt(1.0*n)<1e5){
825 /*
826 if (ii=2){info = DSDPSetY0(dsdp,1,-10000*aa);} else {info = DSDPSetY0(dsdp,m,-10000*aa);}
827 for (i=0;i<n;i++){info = DSDPSetY0(dsdp,i+ii,-100*sqrt(1.0*n)*cnorm);}
828 */
829 /* info=DSDPSetR0(dsdp,cnorm); info=DSDPSetZBar(dsdp,n*n*n*ddmax*cnorm); */ info=DSDPSetPotentialParameter(dsdp,5.0); info=DSDPReuseMatrix(dsdp,0);
830 }
831 }
832 return 0;
833 }
834
835 #undef __FUNCT__
836 #define __FUNCT__ "TCheckArgs0"
TCheckArgs0(DSDP dsdp,SDPCone sdpcone,int m,int nargs,char * runargs[])837 static int TCheckArgs0(DSDP dsdp,SDPCone sdpcone, int m,int nargs,char *runargs[]){
838
839 int kk,info,iloginfo=0;
840 for (kk=1; kk<nargs-1; kk++){
841 if (0){
842 } else if (strncmp(runargs[kk],"-dloginfo",8)==0){
843 iloginfo=atoi(runargs[kk+1]);
844 }
845 }
846 info=DSDPLogInfoAllow(iloginfo,0);
847 return 0;
848 }
849
850 #undef __FUNCT__
851 #define __FUNCT__ "TCheckArgs"
TCheckArgs(DSDP dsdp,SDPCone sdpcone,int m,int nargs,char * runargs[])852 static int TCheckArgs(DSDP dsdp,SDPCone sdpcone, int m,int nargs,char *runargs[]){
853
854 int i,kk, info;
855 int printlevel=10;
856 double *yy0;
857
858 for (kk=1; kk<nargs-1; kk++){
859 if (strncmp(runargs[kk],"-y0",7)==0){
860 DSDPCALLOC2(&yy0,double,m,&info);
861 for (i=0;i<m;i++) yy0[i]=0.0;
862 info = ReadInitialPoint(runargs[kk+1],m,yy0);
863 for (i=0;i<m;i++){info = DSDPSetY0(dsdp,i+1,yy0[i]);}
864 DSDPFREE(&yy0,&info);
865 } else if (strncmp(runargs[kk],"-params",6)==0){
866 info=DSDPReadOptions(dsdp,runargs[kk+1]);
867 } else if (strncmp(runargs[kk],"-print",6)==0){
868 printlevel=atoi(runargs[kk+1]);
869 }
870 }
871
872 info=DSDPSetOptions(dsdp,runargs,nargs);
873 /* if (lpcone){info=LPConeScaleBarrier(lpcone,lpb); */
874 if (rank==0){info=DSDPSetStandardMonitor(dsdp,printlevel);}
875 if (rank==0){info=DSDPSetFileMonitor(dsdp,printlevel);}
876 return(0);
877 }
878
879 #undef __FUNCT__
880 #define __FUNCT__ "ReadInitialPoint"
ReadInitialPoint(char * filename,int m,double yy0[])881 static int ReadInitialPoint(char* filename, int m, double yy0[])
882 {
883 FILE *fp;
884 int i,count;
885
886 fp = fopen(filename,"r");
887 if(!fp)
888 { printf("\n Cannot open file containing initial dual solution %s",filename);
889 return 0;
890 }
891
892 for(count=0,i=0;(i < m) &&(!feof(fp));i++){
893 if(fscanf(fp,"%lf",&yy0[i] ) ){ count++; }
894 }
895
896 if (count < m){
897 printf("Warning: reading initial y vector: \n");
898 printf(" Expecting %d entries but read only %d entries \n",m,count);
899 }
900 fclose(fp);
901 return 0;
902 } /* ReadUserY */
903