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";
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*);
39 int *block,*constraint,*matind;
46 int m;
int n;
int nblocks;
47 int lpn,lpspot,lpblock,lpnnz;
51 double fixedvard,xout;
54static int ReadSDPA2(
char*,DSDPData*);
55static int GetMarkers(
int,
int,
int*,
int*,
int*);
56static int ComputeY0(
DSDP,DSDPData);
58int ReadSDPAFile(
int argc,
char *argv[]);
60#define CHKDATA(a,b,c) { if (c){ printf("Possible problem in variable %d, block %d. \n",a+1,b+1);} }
64#define __FUNCT__ "main"
65int main(
int argc,
char *argv[]){
67 info=ReadSDPAFile(argc,argv);
72#define __FUNCT__ "ReadSDPAFile"
80int ReadSDPAFile(
int argc,
char *argv[]){
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;
99 int *ittt,sspot,ione=1;
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++){
106 if (strncmp(argv[i],
"-benchmark",8)==0){
107 strncpy(thisline,argv[i+1],90); fp1=fopen(thisline,
"r");runbenchmark=1; justone=0;
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;};
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");
123 while ((runbenchmark && !feof(fp1)) || justone==1){
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);
130 strncpy(filename,argv[1],90);
131 strncpy(problemname,argv[1],90);
134 if (rank==0 && fileout){
135 dsdpoutputfile=fopen(outputfile,
"a");
136 fprintf(dsdpoutputfile,
"%s\n",problemname);
137 }
else {dsdpoutputfile=0;fileout=0;}
140 info=ReadSDPA2(filename, &dddd); m=dddd.m;
141 if (info){ printf(
"Problem reading SDPA file\n");
return 1;}
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);
151 for (i=0; i<argc-1; i++){
152 if (strncmp(argv[i],
"-dloginfo",8)==0){ info=DSDPLogInfoAllow(atoi(argv[i+1]),0);};
155 info = DSDPCreate(dddd.m,&dsdp);
156 info = DSDPCreateSDPCone(dsdp,dddd.nblocks,&sdpcone);
158 for (i=0;i<m;i++){info = DSDPSetDualObjective(dsdp,i+1,dddd.dobj[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);
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);
171 spot=0;ijnnz=0;np=0;sdpnmax=1;sdpn=0;stat1=1;
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);
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);
187 }
else if ( ijnnz==0 ){
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 ){
192 info=SDPConeSetADenseVecMat(sdpcone,j,i,n,1.0,dddd.nnz+spot,ijnnz);CHKDATA(i,j,info);
194 info=SDPConeSetASparseVecMat(sdpcone,j,i,n,1.0,0,dddd.matind+spot,dddd.nnz+spot,ijnnz);CHKDATA(i,j,info);
196 if (0==1){ info=SDPConeViewDataMatrix(sdpcone,j,i);}
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];
207 DSDPCALLOC2(&ittt,
int,(m+2),&info);
208 for (i=0;i<=m;i++){ittt[i]=0;}
210 info=GetMarkers(j+1,i,dddd.block+spot,dddd.constraint+spot,&ijnnz);
211 ittt[i+1]=ijnnz; spot+=ijnnz;
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);}
222 info=DSDPCreateBCone(dsdp, &bcone);
223 info=BConeAllocateBounds(bcone,2*m);
225 info=BConeSetUpperBound(bcone,i+1,10);
228 info=BConeSetLowerBound(bcone,i+1,-10);
233 if (printsummary && rank==0){printf(
"DSDP Set Data: %4.3e seconds\n",t3-t2);}
236 if (np<100 && its==0) its=1;
239 if (m<2000 && its>10) its=10;
242 info=DSDPReuseMatrix(dsdp,its);
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);
252 info=DSDPGetDataNorms(dsdp, dnorm);
254 info=DSDPSetR0(dsdp,np);
255 info=DSDPSetGapTolerance(dsdp,1e-3);
256 info=DSDPSetYBounds(dsdp,-1.0,1.0);
259 info = TCheckArgs0(dsdp,sdpcone,dddd.m,argc,argv);
260 info = TCheckArgs(dsdp,sdpcone,dddd.m,argc,argv);
262 info = DSDPSetup(dsdp);
if (info){ printf(
"\nProblem Setting problem. Likely insufficient memory\n");
return 1;}
263 if (0==1){info=SDPConeCheckData(sdpcone);}
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);
277 if (0==1){info=DSDPPrintData(dsdp,sdpcone,lpcone);}
279 info = DSDPSolve(dsdp);
280 if (info){ printf(
"\nNumerical errors encountered in DSDPSolve(). \n");}
282 info=DSDPStopReason(dsdp,&reason);
284 info=DSDPComputeX(dsdp);DSDPCHKERR(info);
286 info=DSDPStopReason(dsdp,&reason);
287 info=DSDPGetSolutionType(dsdp,&pdfeasible);
291 info=DSDPGetDObjective(dsdp,&ddobj);
292 info=DSDPGetPObjective(dsdp,&ppobj);
293 info=DSDPGetFinalErrors(dsdp,derr);
294 info=DSDPGetIts(dsdp,&its);
297 if (printsummary && rank==0){
300 printf(
"DSDP Converged. \n");
303 printf(
"DSDP Terminated Because Dual Objective Exceeded its Bound\n");
306 printf(
"DSDP Terminated Due to Small Steps\n");
309 printf(
"DSDP Terminated Due Maximum Number of Iterations\n");
312 printf(
"DSDP Terminated Due to Infeasible Starting Point\n");
315 printf(
"DSDP Terminated Due to Indefinite Schur Complement\n");
318 printf(
"DSDP Finished\n");
323 printf(
"DSDP Dual Unbounded, Primal Infeasible\n");
325 printf(
"DSDP Primal Unbounded, Dual Infeasible\n");
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);
339 fp2=fopen(tablename,
"a");
341 fprintf(fp2,
" %-18s & %4d & %4d & infeasible & unbounded & %4.0e & %c & %3d & %6.2f \\\\\n",problemname,m,np,derr[0],success,its,t5-t3);
343 fprintf(fp2,
" %-18s & %4d & %4d & unbounded & infeasible & %4.0e & %c & %3d & %6.2f \\\\\n",problemname,m,np,derr[0],success,its,t5-t3);
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);
350 if (printsummary && rank==0){
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);
366 if (printsummary){ DSDPEventLogSummary();}
367 printf(
"--- DSDP Finished ---\n\n");
372 fout=fopen(savefile,
"w");
374 info= DSDPPrintSolution(fout,dsdp,sdpcone,lpcone);
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));
385 for (i=0; i<argc; i++){
386 if (rank==0 && strncmp(argv[i],
"-view",5)==0){DSDPView(dsdp);}
389 info = DSDPDestroy(dsdp);
390 DSDPFREE(&dddd.matind,&info);
391 DSDPFREE(&dddd.nnz,&info);
393 if (fileout){fclose(dsdpoutputfile);}
394 if (0){ DSDPMemoryLog();}
397 if (runbenchmark){ fclose(fp1);}
404#define BUFFERSIZ 4000
405#define BUFFERSIZ 4000
407#define __FUNCT__ "ReadSDPA2"
408static int ReadSDPA2(
char *filename, DSDPData*ddd){
410 char ctmp,refline[BUFFERSIZ]=
"*",thisline[BUFFERSIZ]=
"*";
411 int info,tline,line=0;
415 int i1,nblk,nmat,col,row;
420 fp=fopen(filename,
"r");
422 printf(
"Cannot open file %s !",filename);
return(1);
426 while(!feof(fp) && (thisline[0] ==
'*' || thisline[0] ==
'"') ){
427 fgets(thisline,BUFFERSIZ,fp); line++;
430 if (sscanf(thisline,
"%d",&m)<1){
431 printf(
"Error: line %d. Number of constraints not given.\n",line);
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);
440 ddd->lpn=0;ddd->lpspot=0;ddd->lpblock=0;ddd->cnorm=0;
442 DSDPCALLOC2(&ddd->sformat,
char, (nblocks+1),&info);
443 DSDPCALLOC2(&ddd->blocksizes,
int, (nblocks+1),&info);
444 DSDPCALLOC2(&ddd->conetypes,
char, (nblocks+1),&info );
446 for (i=0;i<nblocks; i++){
447 if (fscanf(fp,
"{")==1 || fscanf(fp,
"(")==1 || fscanf(fp,
",")==1 ){
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';}
456 else{ printf(
"Error block sizes, line %d",line);
return(1);}
458 if (ddd->blocksizes[nblocks-1]==0) nblocks--;
459 fgets(thisline,BUFFERSIZ,fp);
462 DSDPCALLOC2(&ddd->y0,
double,m,&info);
463 DSDPCALLOC2(&ddd->dobj,
double,m,&info);
466 if (fscanf(fp,
",")==1){
470 while (fscanf(fp,
"%lg",&val)!=1){
471 fscanf(fp,
"%c",&ctmp);
473 printf(
"Constraints: %d, Blocks: %d\n",m,nblocks);
474 printf(
"Error reading objective, line %d, i=%d \n",line,i);
return 1;
480 fgets(thisline,BUFFERSIZ,fp);
483 fseek(fp,0,SEEK_SET);
485 for (i=0;i<tline;i++){ctmp=
'*';
while (ctmp!=
'\n'){ fscanf(fp,
"%c",&ctmp);} line++;}
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){
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);
503 printf(
"Data Error in line: %d. Row %d or col %d <= 0 \n%s",line,row+1,col+1,thisline);
506 if (nmat>m || nmat<0){
507 printf(
"Data Error in line: %d. Is Var 0 <= %d <= %d \n%s",line,nmat,m,thisline);
510 if (nblk>nblocks || nblk<0){
511 printf(
"Data Error in line: %d. Is Block 0 <= %d <= %d \n%s",line,nmat,m,thisline);
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){
520 printf(
"Problem Reading SDPA file at line %d: %s\n",line, thisline);
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);
532 fseek(fp,0,SEEK_SET);
534 for (i=0;i<tline;i++){ctmp=
'*';
while (ctmp!=
'\n') fscanf(fp,
"%c",&ctmp); line++;}
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){
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");
550 if (nargs==5 && val!=0.0){
552 printf(
"Warning: Line: %d Row < Column. %s \n",line,thisline);
555 n=ddd->blocksizes[nblk-1];
556 if (nmat==0) {val=-val;}
557 if (ddd->conetypes[nblk-1]==
'S'){
559 ddd->matind[k]=row*(row+1)/2 + col;
560 if (ddd->sformat[nblk-1]==
'U'){ddd->matind[k]=row*n + col;}
565 ddd->constraint[k]=nmat;
568 }
else if (nargs==5 && val==0.0){
569 }
else if (nargs==0){
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");
578 ddd->block[k]=nblocks+1; ddd->constraint[k]=m+2;
579 ddd->matind[k]=10000000; ddd->nnz[k]=0.0;
581 qusort(ddd->block,ddd->constraint,ddd->matind,ddd->nnz,0,nonzero-1);
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];
597 ddd->constraint[nonzero-1]=bigint;ddd->nnz[nonzero-1]=0;
602 ddd->fixedvari=0;ddd->fixedvard=0;
605 if (ddd->blocksizes[ddd->lpblock]==2){
607 while (ddd->block[i]<=ddd->lpblock && i<nonzero){ i++;} spot=i;
608 while (ddd->block[i]==ddd->lpblock+1 && i<nonzero){ i++;k++;}
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;
624 ddd->totalnonzeros=nonzero;
625 for (ddd->n=0,i=0;i<nblocks;i++) ddd->n += ddd->blocksizes[i];
627 ddd->nblocks=nblocks;
634#define __FUNCT__ "Parseline"
635static int Parseline(
char thisline[],
int *nmat,
int *nblk,
int *row,
636 int *col,
double *value,
int *nargs){
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;
652static int partition(
int list1[],
int list2[],
int list3[],
double list5[],
int lstart,
int lend){
654 int pivot1=list1[k],pivot2=list2[k],pivot3=list3[k];
655 double pivot5 = list5[k];
656 int bottom = lstart-1, top = lend;
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];
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];
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]) ){
717 if (ordered && lend-lstart>2){
725static int qusort(
int list1[],
int list2[],
int list3[],
double list5[],
int lstart,
int 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);
740#define __FUNCT__ "GetMarkers"
741static int GetMarkers(
int block,
int constraint,
int blockn[],
int constraintn[],
744 while (blockn[i]==block && constraintn[i]==constraint){ i++;}
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++;}
762#define __FUNCT__ "CheckForConstantMat"
763static int CheckForConstantMat(
double v[],
int nnz,
int n){
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;}
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;
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]);}
784 for (i=1;i<=m;i++,spot+=ijnnz){
785 info=GetMarkers(1,i,dddd.block+spot,dddd.constraint+spot,&ijnnz);
788 dd=dddd.nnz[spot]; ddiag=1;bb=dddd.dobj[i-1];
789 for (ii=0;ii<n;ii++){
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;}
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);
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;
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; }
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);
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;
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; }
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;}
842 if (diag && cnorm*sqrt(1.0*n)<1e5){
847 info=DSDPSetPotentialParameter(dsdp,5.0); info=DSDPReuseMatrix(dsdp,0);
854#define __FUNCT__ "TCheckArgs0"
855static int TCheckArgs0(
DSDP dsdp,
SDPCone sdpcone,
int m,
int nargs,
char *runargs[]){
857 int kk,info,iloginfo=0;
858 for (kk=1; kk<nargs-1; kk++){
860 }
else if (strncmp(runargs[kk],
"-dloginfo",8)==0){
861 iloginfo=atoi(runargs[kk+1]);
864 info=DSDPLogInfoAllow(iloginfo,0);
869#define __FUNCT__ "TCheckArgs"
870static int TCheckArgs(
DSDP dsdp,
SDPCone sdpcone,
int m,
int nargs,
char *runargs[]){
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]);
890 info=DSDPSetOptions(dsdp,runargs,nargs);
892 if (rank==0){info=DSDPSetStandardMonitor(dsdp,printlevel);}
893 if (rank==0){info=DSDPSetFileMonitor(dsdp,printlevel);}
898#define __FUNCT__ "ReadInitialPoint"
899static int ReadInitialPoint(
char* filename,
int m,
double yy0[])
904 fp = fopen(filename,
"r");
906 { printf(
"\n Cannot open file containing initial dual solution %s",filename);
910 for(count=0,i=0;(i < m) &&(!feof(fp));i++){
911 if(fscanf(fp,
"%lf",&yy0[i] ) ){ count++; }
915 printf(
"Warning: reading initial y vector: \n");
916 printf(
" Expecting %d entries but read only %d entries \n",m,count);
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).
struct LPCone_C * LPCone
The LPCone object points to blocks of data that specify linear scalar inequality constraints.
DSDPTerminationReason
There are many reasons to terminate the solver.
@ DSDP_INDEFINITE_SCHUR_MATRIX
DSDPSolutionType
Formulations (P) and (D) can be feasible and bounded, feasible and unbounded, or infeasible.
Internal structures for the DSDP solver.
Internal structure for semidefinite cone.