1 /*
2 * This is an easy to call version of the sdp routine. It takes as
3 * input a problem (n,k,C,a,constraints,constant_offset), and an
4 * initial solution (X,y,Z), allocates working storage, and calls sdp()
5 * to solve the problem. The solution is returned in X,y,Z,pobj,dobj, and
6 * the return code from sdp is returned as the return value from easy_sdp.
7 *
8 */
9
10 #include <stdlib.h>
11 #include <stdio.h>
12 #include <math.h>
13 #include "declarations.h"
14
easy_sdp(n,k,C,a,constraints,constant_offset,pX,py,pZ,ppobj,pdobj)15 int easy_sdp(n,k,C,a,constraints,constant_offset,pX,py,pZ,ppobj,pdobj)
16 int n;
17 int k;
18 struct blockmatrix C;
19 double *a;
20 struct constraintmatrix *constraints;
21 double constant_offset;
22 struct blockmatrix *pX;
23 double **py;
24 struct blockmatrix *pZ;
25 double *ppobj;
26 double *pdobj;
27 {
28 int ret;
29 struct constraintmatrix fill;
30 struct paramstruc params;
31 struct blockmatrix work1;
32 struct blockmatrix work2;
33 struct blockmatrix work3;
34 struct blockmatrix bestx;
35 struct blockmatrix bestz;
36 struct blockmatrix Zi;
37 struct blockmatrix dZ;
38 struct blockmatrix dX;
39 struct blockmatrix cholxinv;
40 struct blockmatrix cholzinv;
41 double *workvec1;
42 double *workvec2;
43 double *workvec3;
44 double *workvec4;
45 double *workvec5;
46 double *workvec6;
47 double *workvec7;
48 double *workvec8;
49 double *diagO;
50 double *Fp;
51 double *O;
52 double *dy;
53 double *dy1;
54 double *rhs;
55 double *besty;
56 int printlevel;
57 int ldam;
58 struct sparseblock **byblocks;
59 struct sparseblock *ptr;
60 struct sparseblock *oldptr;
61 int i;
62 int j;
63 struct sparseblock *p;
64 struct sparseblock *q;
65 struct sparseblock *prev=NULL;
66 double gap;
67 int nnz;
68 int denseblocks;
69 int numblocks;
70
71 /*
72 * Initialize the parameters.
73 */
74
75 initparams(¶ms,&printlevel);
76
77
78 /*
79 * Allocate working storage
80 */
81
82 alloc_mat(C,&work1);
83 alloc_mat(C,&work2);
84 alloc_mat(C,&work3);
85 alloc_mat_packed(C,&bestx);
86 alloc_mat_packed(C,&bestz);
87 alloc_mat_packed(C,&cholxinv);
88 alloc_mat_packed(C,&cholzinv);
89
90 besty=(double *)malloc(sizeof(double)*(k+1));
91 if (besty == NULL)
92 {
93 printf("Storage Allocation Failed!\n");
94 exit(205);
95 };
96
97 if (n > k)
98 {
99 workvec1=(double *)malloc(sizeof(double)*(n+1));
100 }
101 else
102 {
103 workvec1=(double *)malloc(sizeof(double)*(k+1));
104 };
105 if (workvec1 == NULL)
106 {
107 printf("Storage Allocation Failed!\n");
108 exit(205);
109 };
110
111 if (n > k)
112 {
113 workvec2=(double *)malloc(sizeof(double)*(n+1));
114 }
115 else
116 {
117 workvec2=(double *)malloc(sizeof(double)*(k+1));
118 };
119 if (workvec2 == NULL)
120 {
121 printf("Storage Allocation Failed!\n");
122 exit(205);
123 };
124
125 if (n > k)
126 {
127 workvec3=(double *)malloc(sizeof(double)*(n+1));
128 }
129 else
130 {
131 workvec3=(double *)malloc(sizeof(double)*(k+1));
132 };
133 if (workvec3 == NULL)
134 {
135 printf("Storage Allocation Failed!\n");
136 exit(205);
137 };
138
139 if (n > k)
140 {
141 workvec4=(double *)malloc(sizeof(double)*(n+1));
142 }
143 else
144 {
145 workvec4=(double *)malloc(sizeof(double)*(k+1));
146 };
147 if (workvec4 == NULL)
148 {
149 printf("Storage Allocation Failed!\n");
150 exit(205);
151 };
152
153 if (n > k)
154 {
155 workvec5=(double *)malloc(sizeof(double)*(n+1));
156 }
157 else
158 {
159 workvec5=(double *)malloc(sizeof(double)*(k+1));
160 };
161 if (workvec5 == NULL)
162 {
163 printf("Storage Allocation Failed!\n");
164 exit(205);
165 };
166
167 if (n > k)
168 {
169 workvec6=(double *)malloc(sizeof(double)*(n+1));
170 }
171 else
172 {
173 workvec6=(double *)malloc(sizeof(double)*(k+1));
174 };
175 if (workvec6 == NULL)
176 {
177 printf("Storage Allocation Failed!\n");
178 exit(205);
179 };
180
181 if (n > k)
182 {
183 workvec7=(double *)malloc(sizeof(double)*(n+1));
184 }
185 else
186 {
187 workvec7=(double *)malloc(sizeof(double)*(k+1));
188 };
189 if (workvec7 == NULL)
190 {
191 printf("Storage Allocation Failed!\n");
192 exit(205);
193 };
194
195 if (n > k)
196 {
197 workvec8=(double *)malloc(sizeof(double)*(n+1));
198 }
199 else
200 {
201 workvec8=(double *)malloc(sizeof(double)*(k+1));
202 };
203 if (workvec8 == NULL)
204 {
205 printf("Storage Allocation Failed!\n");
206 exit(205);
207 };
208
209
210 if (n > k)
211 {
212 diagO=(double *)malloc(sizeof(double)*(n+1));
213 }
214 else
215 {
216 diagO=(double *)malloc(sizeof(double)*(k+1));
217 };
218 if (diagO == NULL)
219 {
220 printf("Storage Allocation Failed!\n");
221 exit(205);
222 };
223
224
225
226 rhs=malloc(sizeof(double)*(k+1));
227 if (rhs == NULL)
228 {
229 printf("Storage Allocation Failed!\n");
230 exit(205);
231 };
232
233 dy=malloc(sizeof(double)*(k+1));
234 if (dy == NULL)
235 {
236 printf("Storage Allocation Failed!\n");
237 exit(205);
238 };
239
240 dy1=malloc(sizeof(double)*(k+1));
241 if (dy1 == NULL)
242 {
243 printf("Storage Allocation Failed!\n");
244 exit(205);
245 };
246
247 Fp=malloc(sizeof(double)*(k+1));
248 if (Fp == NULL)
249 {
250 printf("Storage Allocation Failed!\n");
251 exit(205);
252 };
253
254 /*
255 * Check to make sure that there is at least one constraint.
256 */
257
258 if (k<= 0)
259 {
260 printf("Problem must have at least one constraint.\n");
261 exit(206);
262 };
263
264 /*
265 * Work out the leading dimension for the array. Note that we may not
266 * want to use k itself, for cache issues.
267 */
268 if ((k % 2) == 0)
269 ldam=k+1;
270 else
271 ldam=k;
272
273 O=malloc(sizeof(double)*ldam*ldam);
274 if (O == NULL)
275 {
276 printf("Storage Allocation Failed!\n");
277 exit(205);
278 };
279
280 alloc_mat(C,&Zi);
281 alloc_mat(C,&dZ);
282 alloc_mat(C,&dX);
283
284 /*
285 * Fill in lots of details in the constraints data structure that haven't
286 * necessarily been done before now.
287 */
288
289 /*
290 * Next, setup issparse and NULL out all nextbyblock pointers.
291 */
292
293 for (i=1; i<=k; i++)
294 {
295 p=constraints[i].blocks;
296 while (p != NULL)
297 {
298 /*
299 * First, set issparse. The 0.125 parameter was hand tuned.
300 */
301
302 if (((1.0*k*k*p->numentries*p->numentries) > 0.125*(1.0*k*p->blocksize*p->blocksize*p->blocksize)) && ((p->numentries) > 5))
303 {
304 p->issparse=0;
305 }
306 else
307 {
308 p->issparse=1;
309 };
310
311 if (C.blocks[p->blocknum].blockcategory == DIAG)
312 p->issparse=1;
313
314 /*
315 * Setup the cross links.
316 */
317
318 p->nextbyblock=NULL;
319 p=p->next;
320 };
321 };
322
323 /*
324 * Now, cross link.
325 */
326 for (i=1; i<=k; i++)
327 {
328 p=constraints[i].blocks;
329 while (p != NULL)
330 {
331 if (p->nextbyblock == NULL)
332 {
333 /*
334 * link in the remaining blocks.
335 */
336 for (j=i+1; j<=k; j++)
337 {
338 q=constraints[j].blocks;
339
340 while (q != NULL)
341 {
342 if (q->blocknum == p->blocknum)
343 {
344 if (p->nextbyblock == NULL)
345 {
346 p->nextbyblock=q;
347 q->nextbyblock=NULL;
348 prev=q;
349 }
350 else
351 {
352 prev->nextbyblock=q;
353 q->nextbyblock=NULL;
354 prev=q;
355 };
356 break;
357 };
358 q=q->next;
359 };
360 };
361 };
362 p=p->next;
363 };
364 };
365
366 /*
367 * If necessary, print out information on sparsity of blocks.
368 */
369
370 if (printlevel >= 1)
371 {
372 numblocks=0;
373 denseblocks=0;
374 for (i=1; i<=k; i++)
375 {
376 p=constraints[i].blocks;
377 while (p != NULL)
378 {
379 if (printlevel >= 4)
380 printf("Sparsity of constraint blocks: %d,%d,%d,%d \n",i,p->blocknum,p->issparse,p->numentries);
381
382 if (p->issparse == 0)
383 denseblocks=denseblocks+1;
384
385 numblocks=numblocks+1;
386 p=p->next;
387 };
388 };
389
390 if (printlevel >= 2)
391 printf("Percentage of dense constraint blocks is %f\n",
392 (100.0*denseblocks)/(1.0*numblocks));
393 };
394
395 /*
396 * Allocate space for byblocks pointers.
397 */
398
399 byblocks=(struct sparseblock **)malloc((C.nblocks+1)*sizeof(struct sparseblock *));
400 if (byblocks == NULL)
401 {
402 printf("Storage Allocation Failed!\n");
403 exit(205);
404 };
405
406 for (i=1; i<=C.nblocks; i++)
407 byblocks[i]=NULL;
408
409 /*
410 * Fill in byblocks pointers.
411 */
412
413 for (i=1; i<=k; i++)
414 {
415 ptr=constraints[i].blocks;
416 while (ptr != NULL)
417 {
418 if (byblocks[ptr->blocknum]==NULL)
419 {
420 byblocks[ptr->blocknum]=ptr;
421 };
422 ptr=ptr->next;
423 };
424 };
425
426 /*
427 * Compute "fill". This data structure tells us which elements in the
428 * block diagonal matrix have corresponding elements in one of the
429 * constraints, and which constraint this element first appears in.
430 *
431 */
432
433 makefill(k,C,constraints,&fill,work1,printlevel);
434
435 /*
436 * Compute the nonzero structure of O. This is only used for debugging
437 * at this point.
438 */
439
440 if (printlevel >= 2)
441 {
442 nnz=structnnz(n,k,C,constraints);
443 printf("Structural density of O %d, %e \n",nnz,nnz*1.0/(k*k*1.0));
444 };
445
446
447 /*
448 * Sort entries in the constraints.
449 */
450
451 sort_entries(k,C,constraints);
452
453 /*
454 * Check the symmetry of C.
455 */
456
457 checkc(n,C,printlevel);
458
459 /*
460 * Check constraints.
461 */
462
463 checkconstraints(n,k,C,constraints,printlevel);
464
465 /*
466 * Now, call sdp().
467 */
468
469 ret=sdp(n,k,C,a,constant_offset,constraints,byblocks,fill,*pX,*py,*pZ,
470 cholxinv,cholzinv,ppobj,pdobj,work1,work2,work3,workvec1,
471 workvec2,workvec3,workvec4,workvec5,workvec6,workvec7,workvec8,
472 diagO,bestx,besty,bestz,Zi,O,rhs,dZ,dX,dy,dy1,Fp,
473 printlevel,params);
474
475 if (printlevel >= 1)
476 {
477 if (ret==0)
478 printf("Success: SDP solved\n");
479 if (ret==1)
480 printf("Success: SDP is primal infeasible\n");
481 if (ret==2)
482 printf("Success: SDP is dual infeasible\n");
483 if (ret==3)
484 printf("Partial Success: SDP solved with reduced accuracy\n");
485 if (ret >= 4)
486 printf("Failure: return code is %d \n",ret);
487
488 if (ret==1)
489 {
490 op_at(k,*py,constraints,work1);
491 addscaledmat(work1,-1.0,*pZ,work1);
492 printf("Certificate of primal infeasibility: a'*y=%.5e, ||A'(y)-Z||=%.5e\n",-1.0,Fnorm(work1));
493 };
494
495 if (ret==2)
496 {
497 op_a(k,constraints,*pX,workvec1);
498 printf("Certificate of dual infeasibility: tr(CX)=%.5e, ||A(X)||=%.5e\n",trace_prod(C,*pX),norm2(k,workvec1+1));
499 };
500
501 if ((ret==0) || (ret>=3))
502 {
503 if (printlevel >= 3)
504 {
505 printf("XZ Gap: %.7e \n",trace_prod(*pZ,*pX));
506 gap=*pdobj-*ppobj;
507 printf("Real Gap: %.7e \n",gap);
508 };
509
510 if (printlevel >= 1)
511 {
512 gap=*pdobj-*ppobj;
513
514 printf("Primal objective value: %.7e \n",*ppobj);
515 printf("Dual objective value: %.7e \n",*pdobj);
516 printf("Relative primal infeasibility: %.2e \n",
517 pinfeas(k,constraints,*pX,a,workvec1));
518 printf("Relative dual infeasibility: %.2e \n",
519 dinfeas(k,C,constraints,*py,*pZ,work1));
520 printf("Real Relative Gap: %.2e \n",gap/(1+fabs(*pdobj)+fabs(*ppobj)));
521 printf("XZ Relative Gap: %.2e \n",trace_prod(*pZ,*pX)/(1+fabs(*pdobj)+fabs(*ppobj)));
522
523 printf("DIMACS error measures: %.2e %.2e %.2e %.2e %.2e %.2e\n",
524 pinfeas(k,constraints,*pX,a,workvec1)*(1+norm2(k,a+1))/
525 (1+norminf(k,a+1)),
526 0.0,
527 dimacserr3(k,C,constraints,*py,*pZ,work1),
528 0.0,
529 gap/(1+fabs(*pdobj)+fabs(*ppobj)),
530 trace_prod(*pZ,*pX)/(1+fabs(*pdobj)+fabs(*ppobj)));
531 };
532 };
533
534 };
535
536 /*
537 * Now, free up all of the storage.
538 */
539
540 free_mat(work1);
541 free_mat(work2);
542 free_mat(work3);
543 free_mat_packed(bestx);
544 free_mat_packed(bestz);
545 free_mat_packed(cholxinv);
546 free_mat_packed(cholzinv);
547
548 free_mat(Zi);
549 free_mat(dZ);
550 free_mat(dX);
551
552 free(besty);
553 free(workvec1);
554 free(workvec2);
555 free(workvec3);
556 free(workvec4);
557 free(workvec5);
558 free(workvec6);
559 free(workvec7);
560 free(workvec8);
561 free(rhs);
562 free(dy);
563 free(dy1);
564 free(Fp);
565 free(O);
566 free(diagO);
567 free(byblocks);
568
569 /*
570 * Free up the fill data structure.
571 */
572
573 ptr=fill.blocks;
574 while (ptr != NULL)
575 {
576 free(ptr->entries);
577 free(ptr->iindices);
578 free(ptr->jindices);
579 oldptr=ptr;
580 ptr=ptr->next;
581 free(oldptr);
582 };
583
584
585 /*
586 * Finally, free the constraints array.
587 */
588
589 return(ret);
590
591 }
592
structnnz(n,k,C,constraints)593 int structnnz(n,k,C,constraints)
594 int n;
595 int k;
596 struct blockmatrix C;
597 struct constraintmatrix *constraints;
598 {
599 int i,j;
600 int ii,jj;
601 int nnz;
602 struct sparseblock *ptri;
603 struct sparseblock *ptrj;
604
605
606 nnz=0;
607 for (i=1; i<=k; i++)
608 for (j=1; j<=k; j++)
609 {
610 ptri=constraints[i].blocks;
611
612 while (ptri != NULL)
613 {
614
615 ptrj=constraints[j].blocks;
616 while (ptrj != NULL)
617 {
618
619 if (ptri->blocknum == ptrj->blocknum)
620 {
621 if (C.blocks[ptri->blocknum].blockcategory==MATRIX)
622 {
623 nnz++;
624 goto NEXTJ;
625 }
626 else
627 { /* DIAG block */
628 for (ii=1; ii<=ptri->numentries; ii++)
629 for (jj=1; jj<=ptrj->numentries; jj++)
630 {
631 if (ptri->iindices[ii]==ptrj->iindices[jj])
632 {
633 nnz++;
634 goto NEXTJ;
635 };
636 };
637 };
638 };
639 ptrj=ptrj->next;
640 };
641
642 ptri=ptri->next;
643 }; /* end while */
644
645 NEXTJ:;
646 }; /* end nested fors */
647
648 return(nnz);
649 }
650
actnnz(n,lda,A)651 int actnnz(n,lda,A)
652 int n;
653 int lda;
654 double A[];
655 {
656 int i,j;
657 int nnz;
658
659 nnz=0;
660 for (i=1; i<=n; i++)
661 {
662 if (A[ijtok(i,i,lda)] != 0.0)
663 nnz++;
664 for (j=i+1; j<=n; j++)
665 {
666 if (A[ijtok(i,j,lda)] != 0.0)
667 {
668 nnz++;
669 nnz++;
670 };
671 };
672 };
673
674 return(nnz);
675 }
676
bandwidth(n,lda,A)677 int bandwidth(n,lda,A)
678 int n;
679 int lda;
680 double A[];
681 {
682 int i;
683 int j;
684 int bw;
685
686 bw=0;
687 for (j=2; j<=n; j++)
688 {
689 for (i=1; i<=j-1; i++)
690 {
691 if (A[ijtok(i,j,lda)] != 0.0)
692 {
693 if ((j-i) > bw)
694 bw=j-i;
695 break;
696 };
697 };
698 };
699
700 return(bw);
701 }
702
703 /*
704 * Sanity checks for the C matrix data structure. If this fails, we'll just
705 * exit(206) to indicate that there was improper input. Otherwise, we'll
706 * return 0;
707 */
708
checkc(int n,struct blockmatrix C,int printlevel)709 int checkc(int n,struct blockmatrix C,int printlevel)
710 {
711 int i,j,k;
712 int totalsize;
713
714 totalsize=0;
715 for (k=1; k<=C.nblocks; k++)
716 {
717 if (C.blocks[k].blockcategory==DIAG)
718 {
719 if (printlevel > 5)
720 printf("blockcategory=diag\n");
721 };
722 if (C.blocks[k].blockcategory==MATRIX)
723 {
724 if (printlevel > 5)
725 printf("blockcategory=matrix\n");
726 };
727
728 totalsize=totalsize+C.blocks[k].blocksize;
729
730 if (C.blocks[k].blockcategory==MATRIX)
731 {
732 for (i=1; i<=C.blocks[k].blocksize; i++)
733 {
734 for (j=1; j<=C.blocks[k].blocksize; j++)
735 {
736 if (C.blocks[k].data.mat[ijtok(i,j,C.blocks[k].blocksize)] !=
737 C.blocks[k].data.mat[ijtok(j,i,C.blocks[k].blocksize)])
738 {
739 if (printlevel >= 1)
740 printf("C is not symmetric, %d, %d, %d\n",k,i,j);
741 exit(206);
742 }
743 };
744 };
745 };
746 };
747
748 if (totalsize != n)
749 {
750 if (printlevel >= 1)
751 printf("Sum of block sizes does not equal n!\n");
752 exit(206);
753 };
754
755 return(0);
756 }
757
758 /*
759 * Sanity tests on the constraints data structure.
760 */
761
checkconstraints(n,k,C,constraints,printlevel)762 int checkconstraints(n,k,C,constraints,printlevel)
763 int n;
764 int k;
765 struct blockmatrix C;
766 struct constraintmatrix *constraints;
767 int printlevel;
768 {
769 int i,j;
770 struct sparseblock *p;
771
772 for (i=1; i<=k; i++)
773 {
774 p=constraints[i].blocks;
775 if (p==NULL)
776 {
777 if (printlevel >= 1)
778 printf("Constraint %d is empty!\n",i);
779 exit(206);
780 };
781 while (p != NULL)
782 {
783 if (p->constraintnum != i)
784 {
785 if (printlevel >= 1)
786 printf("p->constraintnum != i, i=%d \n",i);
787 exit(206);
788 };
789 if (p->blocksize != C.blocks[p->blocknum].blocksize)
790 {
791 if (printlevel >= 1)
792 printf("p->blocksize is wrong, constraint %d \n",i);
793 exit(206);
794 };
795 if (printlevel > 5)
796 printf("Constraint %d, block %d, entries %d\n",i,p->blocknum,p->numentries);
797 for (j=1; j<=p->numentries; j++)
798 {
799
800 /*
801 * For debugging, print out the constraint entry.
802 */
803
804 if (printlevel >6)
805 printf(" (%d, %d)=%lf\n",p->iindices[j],p->jindices[j],p->entries[j]);
806
807 /*
808 * Make sure that all indices are in range.
809 */
810
811 if (p->iindices[j] > C.blocks[p->blocknum].blocksize)
812 {
813 if (printlevel >= 1)
814 printf("i index is larger than blocksize!\n");
815 exit(206);
816 };
817
818 if (p->jindices[j] > C.blocks[p->blocknum].blocksize)
819 {
820 if (printlevel >= 1)
821 printf("j index is larger than blocksize!\n");
822 exit(206);
823 };
824
825 if (p->iindices[j] < 1)
826 {
827 if (printlevel >= 1)
828 printf("i index is less than 1!\n");
829 exit(206);
830 };
831
832 if (p->jindices[j] < 1)
833 {
834 if (printlevel >= 1)
835 printf("j index is less than 1!\n");
836 exit(206);
837 };
838
839 /*
840 * Make sure that entries are sorted properly.
841 */
842
843 if (p->iindices[j] > p->jindices[j])
844 {
845 if (printlevel >= 1)
846 {
847 printf("i index is greater than j index!\n");
848 printf("constraint=%d\n",i);
849 printf("iindex=%d\n",p->iindices[j]);
850 printf("jindex=%d\n",p->jindices[j]);
851 };
852
853 exit(206);
854 };
855
856
857 /*
858 * Check for any duplicate entries that snuck in.
859 */
860
861 if (j < p->numentries)
862 if ((p->iindices[j]==p->iindices[j+1]) &
863 (p->jindices[j]==p->jindices[j+1]))
864 {
865 if (printlevel >= 1)
866 {
867 printf("Duplicate entry!\n");
868 printf("constraint=%d\n",i);
869 printf("iindex=%d\n",p->iindices[j]);
870 printf("jindex=%d\n",p->jindices[j]);
871 };
872
873 exit(206);
874 };
875
876 if (C.blocks[p->blocknum].blockcategory==DIAG)
877 {
878 if (p->iindices[j] != p->jindices[j])
879 {
880 if (printlevel >= 1)
881 printf("Off diagonal entry in diagonal block!\n");
882 exit(206);
883 };
884 };
885 };
886 p=p->next;
887 };
888 };
889 return(0);
890 }
891