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