19 #ifndef LIMBO_SOLVERS_CSDPEASYSDPAPI_H
20 #define LIMBO_SOLVERS_CSDPEASYSDPAPI_H
32 #include <limbo/thirdparty/Csdp/include/declarations.h>
65 struct constraintmatrix *constraints,
66 double constant_offset,
67 struct blockmatrix *pX,
69 struct blockmatrix *pZ,
73 struct paramstruc
const& params,
78 struct constraintmatrix fill;
80 struct blockmatrix work1;
81 struct blockmatrix work2;
82 struct blockmatrix work3;
83 struct blockmatrix bestx;
84 struct blockmatrix bestz;
85 struct blockmatrix Zi;
86 struct blockmatrix dZ;
87 struct blockmatrix dX;
88 struct blockmatrix cholxinv;
89 struct blockmatrix cholzinv;
107 struct sparseblock **byblocks;
108 struct sparseblock *ptr;
109 struct sparseblock *oldptr;
113 struct sparseblock *p;
114 struct sparseblock *q;
115 struct sparseblock *prev=NULL;
132 alloc_mat_packed(C,&bestx);
133 alloc_mat_packed(C,&bestz);
134 alloc_mat_packed(C,&cholxinv);
135 alloc_mat_packed(C,&cholzinv);
137 besty=(
double *)malloc(
sizeof(
double)*(k+1));
140 printf(
"Storage Allocation Failed!\n");
146 workvec1=(
double *)malloc(
sizeof(
double)*(n+1));
150 workvec1=(
double *)malloc(
sizeof(
double)*(k+1));
152 if (workvec1 == NULL)
154 printf(
"Storage Allocation Failed!\n");
160 workvec2=(
double *)malloc(
sizeof(
double)*(n+1));
164 workvec2=(
double *)malloc(
sizeof(
double)*(k+1));
166 if (workvec2 == NULL)
168 printf(
"Storage Allocation Failed!\n");
174 workvec3=(
double *)malloc(
sizeof(
double)*(n+1));
178 workvec3=(
double *)malloc(
sizeof(
double)*(k+1));
180 if (workvec3 == NULL)
182 printf(
"Storage Allocation Failed!\n");
188 workvec4=(
double *)malloc(
sizeof(
double)*(n+1));
192 workvec4=(
double *)malloc(
sizeof(
double)*(k+1));
194 if (workvec4 == NULL)
196 printf(
"Storage Allocation Failed!\n");
202 workvec5=(
double *)malloc(
sizeof(
double)*(n+1));
206 workvec5=(
double *)malloc(
sizeof(
double)*(k+1));
208 if (workvec5 == NULL)
210 printf(
"Storage Allocation Failed!\n");
216 workvec6=(
double *)malloc(
sizeof(
double)*(n+1));
220 workvec6=(
double *)malloc(
sizeof(
double)*(k+1));
222 if (workvec6 == NULL)
224 printf(
"Storage Allocation Failed!\n");
230 workvec7=(
double *)malloc(
sizeof(
double)*(n+1));
234 workvec7=(
double *)malloc(
sizeof(
double)*(k+1));
236 if (workvec7 == NULL)
238 printf(
"Storage Allocation Failed!\n");
244 workvec8=(
double *)malloc(
sizeof(
double)*(n+1));
248 workvec8=(
double *)malloc(
sizeof(
double)*(k+1));
250 if (workvec8 == NULL)
252 printf(
"Storage Allocation Failed!\n");
259 diagO=(
double *)malloc(
sizeof(
double)*(n+1));
263 diagO=(
double *)malloc(
sizeof(
double)*(k+1));
267 printf(
"Storage Allocation Failed!\n");
273 rhs=(
double*)malloc(
sizeof(
double)*(k+1));
276 printf(
"Storage Allocation Failed!\n");
280 dy=(
double*)malloc(
sizeof(
double)*(k+1));
283 printf(
"Storage Allocation Failed!\n");
287 dy1=(
double*)malloc(
sizeof(
double)*(k+1));
290 printf(
"Storage Allocation Failed!\n");
294 Fp=(
double*)malloc(
sizeof(
double)*(k+1));
297 printf(
"Storage Allocation Failed!\n");
310 O=(
double*)malloc(
sizeof(
double)*ldam*ldam);
313 printf(
"Storage Allocation Failed!\n");
337 p=constraints[i].blocks;
343 if (((p->numentries) > 0.25*(p->blocksize)) && ((p->numentries) > 15))
352 if (C.blocks[p->blocknum].blockcategory == DIAG)
369 p=constraints[i].blocks;
372 if (p->nextbyblock == NULL)
379 for (j=i+1; j<=k; j++)
381 q=constraints[j].blocks;
385 if (q->blocknum == p->blocknum)
387 if (p->nextbyblock == NULL)
417 p=constraints[i].blocks;
420 printf(
"%d,%d,%d,%d \n",i,p->blocknum,p->issparse,p->numentries);
430 byblocks=(
struct sparseblock **)malloc((C.nblocks+1)*
sizeof(
struct sparseblock *));
431 if (byblocks == NULL)
433 printf(
"Storage Allocation Failed!\n");
437 for (i=1; i<=C.nblocks; i++)
445 ptr=constraints[i].blocks;
448 if (byblocks[ptr->blocknum]==NULL)
450 byblocks[ptr->blocknum]=ptr;
463 makefill(k,C,constraints,&fill,work1,printlevel);
469 nnz=structnnz(n,k,C,constraints);
472 printf(
"Structural density of O %d, %e \n",nnz,nnz*1.0/(k*k*1.0));
478 sort_entries(k,C,constraints);
484 ret=sdp(n,k,C,a,constant_offset,constraints,byblocks,fill,*pX,*py,*pZ,
485 cholxinv,cholzinv,ppobj,pdobj,work1,work2,work3,workvec1,
486 workvec2,workvec3,workvec4,workvec5,workvec6,workvec7,workvec8,
487 diagO,bestx,besty,bestz,Zi,O,rhs,dZ,dX,dy,dy1,Fp,
493 printf(
"Success: SDP solved\n");
495 printf(
"Success: SDP is primal infeasible\n");
497 printf(
"Success: SDP is dual infeasible\n");
499 printf(
"Partial Success: SDP solved with reduced accuracy\n");
501 printf(
"Failure: return code is %d \n",ret);
505 op_at(k,*py,constraints,work1);
506 addscaledmat(work1,-1.0,*pZ,work1);
507 printf(
"Certificate of primal infeasibility: a'*y=%.5e, ||A'(y)-Z||=%.5e\n",-1.0,Fnorm(work1));
512 op_a(k,constraints,*pX,workvec1);
513 printf(
"Certificate of dual infeasibility: tr(CX)=%.5e, ||A(X)||=%.5e\n",trace_prod(C,*pX),norm2(k,workvec1+1));
516 if ((ret==0) || (ret>=3))
520 printf(
"XZ Gap: %.7e \n",trace_prod(*pZ,*pX));
522 printf(
"Real Gap: %.7e \n",gap);
529 printf(
"Primal objective value: %.7e \n",*ppobj);
530 printf(
"Dual objective value: %.7e \n",*pdobj);
531 printf(
"Relative primal infeasibility: %.2e \n",
532 pinfeas(k,constraints,*pX,a,workvec1));
533 printf(
"Relative dual infeasibility: %.2e \n",
534 dinfeas(k,C,constraints,*py,*pZ,work1));
535 printf(
"Real Relative Gap: %.2e \n",gap/(1+fabs(*pdobj)+fabs(*ppobj)));
536 printf(
"XZ Relative Gap: %.2e \n",trace_prod(*pZ,*pX)/(1+fabs(*pdobj)+fabs(*ppobj)));
538 printf(
"DIMACS error measures: %.2e %.2e %.2e %.2e %.2e %.2e\n",
539 pinfeas(k,constraints,*pX,a,workvec1)*(1+norm2(k,a+1))/
542 dimacserr3(k,C,constraints,*py,*pZ,work1),
544 gap/(1+fabs(*pdobj)+fabs(*ppobj)),
545 trace_prod(*pZ,*pX)/(1+fabs(*pdobj)+fabs(*ppobj)));
558 free_mat_packed(bestx);
559 free_mat_packed(bestz);
560 free_mat_packed(cholxinv);
561 free_mat_packed(cholzinv);
int easy_sdp_ext(int n, int k, struct blockmatrix C, double *a, struct constraintmatrix *constraints, double constant_offset, struct blockmatrix *pX, double **py, struct blockmatrix *pZ, double *ppobj, double *pdobj, struct paramstruc const ¶ms, int const &printlevel)
API to call Csdp solver.
mathematical utilities such as abs