1 /*
2 Solve a semidefinite programming problem using the algorithm of Helmberg,
3 Rendl, Vanderbei, and Wolkowicz.
4
5 Note: This version of the code uses predictor correct steps.
6
7 Usage:
8
9
10 Return codes:
11
12 0 Success. Problem solved to optimality.
13 1 Success. The problem is primal infeasibile, and we
14 have a certificate.
15 2 Success. The problem is dual infeasible, and we have
16 a certificate.
17
18 3. Partial Success. Didn't reach full accuracy.
19
20 4 Failure: Maximum iterations reached.
21 5 Failure: Stuck at edge of primal feasibility.
22 6 Failure: Stuck at edge of dual feasibility.
23 7 Failure: Lack of progress
24 8 Failure: X, Z, or O was singular.
25 9 Failure: Detected NaN or Inf values.
26 10 Failure: Stopped by QUIT, KILL, TERM, SIGXCPU signal
27
28 Notes on data storage: All "2-d" arrays are stored in Fortran style,
29 column major order. macros in index.h are used to handle indexing into
30 C 1-d vectors. All "1-d" arrays are stored in C vectors, with the first
31 element of the vector (index 0) unused. We always start indexing in
32 both one and two dimensional arrays with 1, ala Fortran. This makes
33 the C code somewhat clearer, but requires us to pass v+1 into a Fortran
34 subroutine instead of just passing v.
35
36 */
37
38 #include <stdio.h>
39 #include <stdlib.h>
40 #include <math.h>
41 #include "declarations.h"
42
43 #ifdef USEGETTIME
44 /*
45 * Stuff for keeping track of time.
46 */
47
48
49 #include <stddef.h> /* definition of NULL */
50 #include <sys/time.h> /* definition of timeval struct and protyping of gettime
51 ofday */
52
53 double opotime=0.0;
54 double factortime=0.0;
55 double totaltime=0.0;
56 double othertime=0.0;
57 double othertime1=0.0;
58 double othertime2=0.0;
59 double othertime3=0.0;
60 struct timeval tp;
61 double t1=0.0;
62 double t2=0.0;
63 double starttime=0.0;
64 double endtime=0.0;
65
66 #endif
67
sdp(n,k,C,a,constant_offset,constraints,byblocks,fill,X,y,Z,cholxinv,cholzinv,pobj,dobj,work1,work2,work3,workvec1,workvec2,workvec3,workvec4,workvec5,workvec6,workvec7,workvec8,diagO,bestx,besty,bestz,Zi,O,rhs,dZ,dX,dy,dy1,Fp,printlevel,parameters)68 int sdp(n,k,C,a,constant_offset,constraints,byblocks,fill,X,y,Z,cholxinv,
69 cholzinv,pobj,dobj,work1,work2,work3,workvec1,workvec2,workvec3,
70 workvec4,workvec5,workvec6,workvec7,workvec8,diagO,bestx,besty,bestz,
71 Zi,O,rhs,dZ,dX,dy,dy1,Fp,printlevel,parameters)
72 int n;
73 int k;
74 struct blockmatrix C;
75 double *a;
76 double constant_offset;
77 struct constraintmatrix *constraints;
78 struct sparseblock **byblocks;
79 struct constraintmatrix fill;
80 struct blockmatrix X;
81 double *y;
82 struct blockmatrix Z;
83 struct blockmatrix cholxinv;
84 struct blockmatrix cholzinv;
85 double *pobj;
86 double *dobj;
87 struct blockmatrix work1;
88 struct blockmatrix work2;
89 struct blockmatrix work3;
90 double *workvec1;
91 double *workvec2;
92 double *workvec3;
93 double *workvec4;
94 double *workvec5;
95 double *workvec6;
96 double *workvec7;
97 double *workvec8;
98 double *diagO;
99 struct blockmatrix bestx;
100 double *besty;
101 struct blockmatrix bestz;
102 struct blockmatrix Zi;
103 double *O;
104 double *rhs;
105 struct blockmatrix dZ;
106 struct blockmatrix dX;
107 double *dy;
108 double *dy1;
109 double *Fp;
110 int printlevel;
111 struct paramstruc parameters;
112 {
113 double gap;
114 double relgap;
115 double mu=1.0e30;
116 double muplus;
117 double muk;
118 double gamma;
119 double alphap,alphap1;
120 double alphad,alphad1;
121 double scale1;
122 double scale2;
123 double nrmx;
124 double relpinfeas=1.0e10;
125 double reldinfeas=1.0e10;
126 double newrelpinfeas;
127 double newreldinfeas;
128 double limrelpinfeas;
129 double bestrelpinfeas;
130 double limreldinfeas;
131 double bestreldinfeas;
132 int iter=0;
133 int m;
134 int ret;
135 int i;
136 int j;
137 int info;
138 int ldam;
139 int retries=0;
140 int retcode=0;
141 double bestmeas;
142 double diagnrm;
143 double diagfact=0.0;
144 double diagadd;
145 double ent;
146
147 /*
148 * Stuff for instrumentation of iterative refinement.
149 */
150 double relerr1;
151 double relerr2=0.0;
152
153 int refinements;
154 double besterr;
155 double relerr;
156
157 /*
158 * Stuff for iterative refinements.
159 */
160
161 int lastimproverefinement;
162 double mindiag;
163
164 /*
165 *
166 */
167
168 int ispfeasprob;
169 int isdfeasprob;
170
171 double normC;
172 double norma;
173
174 /*
175 * Amount by which to perturb
176 */
177
178 double perturbfac;
179
180 /*
181 * Sticky flags for primal and dual feasibility.
182 */
183
184 int probpfeas=0;
185 int probdfeas=0;
186
187 /*
188 *
189 *
190 */
191
192 double pinfeasmeas;
193 double dinfeasmeas;
194
195 /*
196 * Stuff for adjusting the step fraction.
197 */
198
199 double mystepfrac;
200 double minalpha;
201
202 /*
203 *
204 */
205 double newalphap;
206
207 /*
208 * Stuff for keeping track of best solutions.
209 */
210
211 #define BASIZE 100
212 double bestarray[BASIZE+1];
213
214 /*
215 * Used in checking whether the primal-dual affine step gives us a new
216 * best solution.
217 */
218
219 double affgap,affpobj,affdobj,affrelgap,affrelpinfeas,affreldinfeas;
220
221 /*
222 * Indicate the version of CSDP.
223 */
224
225 if (printlevel >= 1)
226 printf("CSDP 6.2.0\n");
227
228 /*
229 * Precompute norms of a and C, so that we don't keep doing this
230 * again and again.
231 *
232 */
233
234 norma=norm2(k,a+1);
235 normC=Fnorm(C);
236 if (parameters.perturbobj>0)
237 /* perturbfac=parameters.perturbobj*1.0e-6*normC/sqrt(1.0*n); */
238 perturbfac=1.0e-6*normC/sqrt(1.0*n);
239 else
240 perturbfac=0.0;
241
242 /*
243 * Determine whether or not this is a feasibility problem.
244 *
245 */
246
247 if (normC==0.0)
248 {
249 if (printlevel >= 1)
250 printf("This is a pure primal feasibility problem.\n");
251 ispfeasprob=1;
252 }
253 else
254 ispfeasprob=0;
255
256 if (norma==0.0)
257 {
258 if (printlevel >= 1)
259 printf("This is a pure dual feasibility problem.\n");
260 isdfeasprob=1;
261 }
262 else
263 isdfeasprob=0;
264
265 /*
266 * For historical reasons, we use m for the total number of constraints.
267 * In this version of the code, k is also the number of constraints.
268 */
269
270 m=k;
271
272 /*
273 * Work out the leading dimension for the array. Note that we may not
274 * want to use k itself, for cache issues.
275 */
276 if ((k % 2) == 0)
277 ldam=k+1;
278 else
279 ldam=k;
280
281 /*
282 * Compute Cholesky factors of X and Z.
283 */
284
285 copy_mat(X,work1);
286 ret=chol(work1);
287
288 /*
289 * If the matrix was singular, then just return with an error.
290 */
291
292 if (ret != 0)
293 {
294 if (printlevel >= 1)
295 printf("X was singular!\n");
296 retcode=8;
297 goto RETURNBEST;
298 };
299
300 chol_inv(work1,work2);
301 store_packed(work2,cholxinv);
302
303 copy_mat(Z,work1);
304 ret=chol(work1);
305
306 /*
307 * If the matrix was singular, then just return with an error.
308 */
309
310 if (ret != 0)
311 {
312 if (printlevel >= 1)
313 printf("Z was singular!\n");
314 retcode=8;
315 goto RETURNBEST;
316 };
317
318 chol_inv(work1,work2);
319 store_packed(work2,cholzinv);
320
321 /*
322 Compute Zi.
323 */
324
325 copy_mat(work2,work1);
326 trans(work1);
327 scale1=1.0;
328 scale2=0.0;
329 mat_mult(scale1,scale2,work2,work1,Zi);
330
331 if (printlevel >= 4)
332 printf("Fnorm of Zi is %e \n",Fnorm(Zi));
333
334 /*
335 Compute primal and dual objective values.
336 */
337
338 *pobj=calc_pobj(C,X,constant_offset);
339 *dobj=calc_dobj(k,a,y,constant_offset);
340 if (parameters.usexzgap==0)
341 {
342 gap=*dobj-*pobj;
343 if (gap < 0.0)
344 gap=0.0;
345 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
346 }
347 else
348 {
349 gap=trace_prod(X,Z);
350 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
351 };
352
353 if (printlevel >= 3)
354 {
355 printf("constant offset is %e \n",constant_offset);
356 printf("Fnorm of X is %e \n",Fnorm(X));
357 printf("Fnorm of C is %e \n",normC);
358 };
359
360 if (printlevel >= 4)
361 {
362 printf("pobj is %e \n",*pobj);
363 printf("dobj is %e \n",*dobj);
364 printf("gap is %e \n",gap);
365 };
366
367 relpinfeas=pinfeas(k,constraints,X,a,workvec1);
368 bestrelpinfeas=relpinfeas;
369 reldinfeas=dinfeas(k,C,constraints,y,Z,work2);
370 bestreldinfeas=reldinfeas;
371
372 if (relpinfeas < parameters.axtol )
373 probpfeas=1;
374
375 if (reldinfeas < parameters.atytol)
376 probdfeas=1;
377
378 /*
379 * Record this solution as the best solution.
380 */
381
382 bestmeas=1.0e100;
383 store_packed(X,bestx);
384 store_packed(Z,bestz);
385 for (i=1; i<=k; i++)
386 besty[i]=y[i];
387
388 /*
389 Initialize the big loop.
390 */
391
392 alphap=0.0;
393 alphad=0.0;
394
395 /*
396 Print out some status information.
397 */
398
399 if (printlevel >= 1)
400 {
401 printf("Iter: %2d Ap: %.2e Pobj: % .7e Ad: %.2e Dobj: % .7e \n",iter,alphap,*pobj,alphad,*dobj);
402 fflush(stdout);
403 };
404
405 while ((relgap > parameters.objtol) || (relpinfeas > parameters.axtol) || (reldinfeas > parameters.atytol))
406 {
407
408
409
410 /*
411 * Call the user exit routine, and let the user stop the process
412 * if he wants to.
413 */
414
415 ret=user_exit(n,k,C,a,*dobj,*pobj,constant_offset,constraints,
416 X,y,Z,parameters);
417 if (ret >= 1)
418 {
419 if (ret==1)
420 {
421 /*
422 * Received SIGTERM, SIGQUIT, SIGXCPU, ...
423 */
424 return(10);
425 }
426 else
427 {
428 /*
429 * User exit was good enough, returning success.
430 */
431 return(0);
432 };
433
434 }
435
436 bestarray[iter % BASIZE]=bestmeas;
437
438 /*
439 * Compute the stepfrac to be used for this iteration.
440 */
441
442 if (iter > 1)
443 {
444 if (alphap > alphad)
445 minalpha=alphad;
446 else
447 minalpha=alphap;
448 }
449 else
450 {
451 minalpha=1.0;
452 };
453
454 mystepfrac=parameters.minstepfrac+minalpha*(parameters.maxstepfrac-parameters.minstepfrac);
455 if (printlevel >= 3)
456 printf("mystepfrac is %e \n",mystepfrac);
457
458 /*
459 * Print out information on primal/dual feasibility.
460 */
461
462 if (printlevel >= 2)
463 {
464 printf("Relative primal infeasibility is %.7e \n",relpinfeas);
465 printf("Relative dual infeasibility: %.7e \n",reldinfeas);
466 printf("Relative duality gap is %.7e \n",relgap);
467 printf("XZ relative duality gap is %.7e \n",trace_prod(X,Z)/(1+fabs(*dobj)));
468 };
469
470 /*
471 * If this is a feasibility problem, and we've got feasibility,
472 * then we're done!
473 */
474
475 if ((ispfeasprob==1) && (relpinfeas < parameters.axtol))
476 {
477 if (printlevel >= 3)
478 printf("Got primal feasibility, so stopping.\n");
479
480 for (i=1; i<=k; i++)
481 y[i]=0.0;
482
483 alphad=parameters.atytol/(sqrt(n*1.0)*200);
484 make_i(work1);
485 if (alphad*trace_prod(X,work1) > parameters.objtol)
486 alphad=0.005*parameters.objtol/trace_prod(X,work1);
487 zero_mat(work2);
488 addscaledmat(work2,alphad,work1,Z);
489
490 relpinfeas=pinfeas(k,constraints,X,a,workvec1);
491 if (relpinfeas < bestrelpinfeas)
492 bestrelpinfeas=relpinfeas;
493 reldinfeas=dinfeas(k,C,constraints,y,Z,work2);
494 if (reldinfeas < bestreldinfeas)
495 bestreldinfeas=reldinfeas;
496
497 *pobj=calc_pobj(C,X,constant_offset);
498 *dobj=calc_dobj(k,a,y,constant_offset);
499
500 if (parameters.usexzgap==0)
501 {
502 gap=*dobj-*pobj;
503 if (gap < 0.0)
504 gap=0.0;
505 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
506 }
507 else
508 {
509 gap=trace_prod(X,Z);
510 if (gap < 0.0)
511 gap=0.0;
512 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
513 };
514
515 bestmeas=relpinfeas/parameters.axtol;
516
517 store_packed(X,bestx);
518 store_packed(Z,bestz);
519 for (i=1; i<=k; i++)
520 besty[i]=y[i];
521
522 retcode=0;
523 goto RETURNBEST;
524
525 };
526
527
528 if ((isdfeasprob==1) && (reldinfeas < parameters.atytol))
529 {
530 if (printlevel >= 3)
531 printf("Got dual feasibility, so stopping.\n");
532
533 make_i(work1);
534 zero_mat(work2);
535 op_a(k,constraints,work1,workvec1);
536 alphap=parameters.atytol/(200.0*norm2(k,workvec1+1));
537 if (alphap*trace_prod(work1,Z) > parameters.objtol)
538 alphap=0.005*parameters.objtol/trace_prod(work1,Z);
539 addscaledmat(work2,alphap,work1,X);
540
541 relpinfeas=pinfeas(k,constraints,X,a,workvec1);
542 if (relpinfeas < bestrelpinfeas)
543 bestrelpinfeas=relpinfeas;
544 reldinfeas=dinfeas(k,C,constraints,y,Z,work2);
545 if (reldinfeas < bestreldinfeas)
546 bestreldinfeas=reldinfeas;
547
548 *pobj=calc_pobj(C,X,constant_offset);
549 *dobj=calc_dobj(k,a,y,constant_offset);
550
551 if (parameters.usexzgap==0)
552 {
553 gap=*dobj-*pobj;
554 if (gap < 0.0)
555 gap=0.0;
556 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
557 }
558 else
559 {
560 gap=trace_prod(X,Z);
561 if (gap < 0.0)
562 gap=0.0;
563 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
564 };
565
566 bestmeas=reldinfeas/parameters.atytol;
567
568 store_packed(X,bestx);
569 store_packed(Z,bestz);
570 for (i=1; i<=k; i++)
571 besty[i]=y[i];
572 if (printlevel >= 3)
573 printf("New best solution, %e \n",bestmeas);
574
575 retcode=0;
576 goto RETURNBEST;
577
578 };
579
580 /*
581 * Check for primal or dual infeasibility.
582 */
583
584
585 op_at(k,y,constraints,work1);
586 addscaledmat(work1,-1.0,Z,work1);
587
588 pinfeasmeas=-(*dobj)/Fnorm(work1);
589 if (printlevel >= 2)
590 printf("-a'*y/||A'(y)-Z|| is %e \n",pinfeasmeas);
591
592 if ((probpfeas==0) &&
593 (pinfeasmeas > parameters.pinftol))
594 {
595 if (printlevel >= 1)
596 printf("Declaring primal infeasibility.\n");
597 for (i=1; i<=k; i++)
598 y[i]=-y[i]/(*dobj);
599 zero_mat(work1);
600 addscaledmat(work1,-1.0/(*dobj),Z,work1);
601 copy_mat(work1,Z);
602 retcode=1;
603 goto RETURNCERT;
604 };
605
606 op_a(k,constraints,X,workvec1);
607 dinfeasmeas=trace_prod(C,X)/norm2(k,workvec1+1);
608 if (printlevel >= 2)
609 printf("<C,X>/||A(X)||=%e\n",dinfeasmeas);
610 if ((probdfeas==0) &&
611 (dinfeasmeas>parameters.dinftol))
612 {
613 if (printlevel >= 1)
614 printf("Declaring dual infeasibility.\n");
615 zero_mat(work1);
616 addscaledmat(work1,1.0/trace_prod(C,X),X,work1);
617 copy_mat(work1,X);
618 retcode=2;
619 goto RETURNCERT;
620 };
621
622 /*
623 * Print out the norm(X) for debugging purposes.
624 */
625
626 if (printlevel >= 3)
627 {
628 nrmx=Fnorm(X);
629 printf("Fnorm of X is %e \n",nrmx);
630 };
631
632
633
634
635 /*
636 Now, compute the system matrix.
637 */
638
639 #ifdef USEGETTIME
640 gettimeofday(&tp,NULL);
641 t1=(double)tp.tv_sec+(1.0e-6)*tp.tv_usec;
642 #endif
643
644 op_o(k,constraints,byblocks,Zi,X,O,work1,work2);
645
646 #ifdef USEGETTIME
647 gettimeofday(&tp,NULL);
648 t2=(double)tp.tv_sec+(1.0e-6)*tp.tv_usec;
649 opotime=opotime+t2-t1;
650 #endif
651
652 /*
653 * Print out the actual density of Z and X.
654 */
655 if ((iter==5) && (printlevel >= 3))
656 {
657 printf("Actual density of O %e\n",actnnz(k,ldam,O)/(1.0*k*k));
658 for (j=1; j<=X.nblocks; j++)
659 {
660 if (X.blocks[j].blockcategory==MATRIX)
661 {
662 printf("density X block %d, %e \n",j,actnnz(X.blocks[j].blocksize,
663 X.blocks[j].blocksize,
664 X.blocks[j].data.mat)/
665 (1.0*X.blocks[j].blocksize*X.blocks[j].blocksize));
666
667 printf("density Z block %d, %e \n",j,actnnz(Z.blocks[j].blocksize,
668 Z.blocks[j].blocksize,
669 Z.blocks[j].data.mat)/
670 (1.0*Z.blocks[j].blocksize*Z.blocks[j].blocksize));
671
672 printf("bandwidth Z block %d, %d/%d \n",j,bandwidth(Z.blocks[j].blocksize,Z.blocks[j].blocksize,Z.blocks[j].data.mat),Z.blocks[j].blocksize);
673
674 };
675 };
676 };
677
678 /*
679 Save a copy of O in the lower diagonal for later use.
680 */
681
682 for (i=1; i<=k-1; i++)
683 for (j=i; j<=k; j++)
684 O[ijtok(j,i, (long int) ldam)]=O[ijtok(i,j, (long int) ldam)];
685
686 for (i=1; i<=k; i++)
687 diagO[i]=O[ijtok(i,i,(long int) ldam)];
688
689 mindiag=1.0e30;
690
691 for (i=1; i<=k; i++)
692 {
693 if (diagO[i] < mindiag)
694 {
695 mindiag=diagO[i];
696 };
697
698 };
699 /*
700 This is where we come if factorization failed or the system
701 was so badly conditioned that we didn't get a usable solution.
702 */
703
704 diagnrm=0.0;
705
706 for (i=1; i<=k; i++)
707 {
708 ent=diagO[i];
709 diagnrm = diagnrm + ent*ent;
710 };
711 diagnrm=sqrt(diagnrm);
712
713
714 RETRYFACTOR:
715
716 /*
717 Now, let's make sure that O isn't singular.
718 */
719
720
721
722
723 diagadd=1.0e-17*diagfact*diagnrm/sqrt(k*1.0);
724
725 while ((diagadd + mindiag) <= 0.0)
726 {
727 retries++;
728 if (diagfact==0.0)
729 {
730 diagfact=0.1;
731 diagadd=1.0e-17*diagfact*diagnrm/sqrt(k*1.0);
732 }
733 else
734 {
735 diagfact=diagfact*10.0;
736 diagadd=diagadd*10.0;
737 };
738 };
739
740 if (printlevel >= 3)
741 printf("diagnrm is %e, adding diagadd %e \n",diagnrm,diagadd);
742
743 for (i=1; i<=k; i++)
744 O[ijtok(i,i, (long int) ldam)] += diagadd;
745
746 /*
747 * Scale the O matrix.
748 */
749
750 for (i=1; i<=k; i++)
751 workvec8[i]=1.0/sqrt(O[ijtok(i,i, (long int) ldam)]);
752
753 for (i=1; i<=k; i++)
754 {
755 if (workvec8[i] > 1.0e30)
756 workvec8[i]=1.0e30;
757 };
758
759 #pragma omp parallel for schedule(dynamic,64) default(none) shared(O,ldam,k,workvec8) private(i,j)
760 for (j=1; j<=k; j++)
761 for (i=1; i<=j; i++)
762 O[ijtok(i,j, (long int) ldam)]=O[ijtok(i,j, (long int) ldam)]*(workvec8[i]*workvec8[j]);
763
764 /*
765 Next, compute the cholesky factorization of the system matrix.
766 */
767
768 #ifdef USEGETTIME
769 gettimeofday(&tp,NULL);
770 t1=(double)tp.tv_sec+(1.0e-6)*tp.tv_usec;
771 #endif
772
773 #ifdef NOUNDERLAPACK
774 #ifdef CAPSLAPACK
775 DPOTRF("U",&m,O,&ldam,&info);
776 #else
777 dpotrf("U",&m,O,&ldam,&info);
778 #endif
779 #else
780 #ifdef CAPSLAPACK
781 DPOTRF_("U",&m,O,&ldam,&info);
782 #else
783 dpotrf_("U",&m,O,&ldam,&info);
784 #endif
785 #endif
786
787 #ifdef USEGETTIME
788 gettimeofday(&tp,NULL);
789 t2=(double)tp.tv_sec+(1.0e-6)*tp.tv_usec;
790 factortime=factortime+t2-t1;
791 #endif
792
793 if (info != 0)
794 {
795 if (printlevel >= 3)
796 printf("Factorization of the system matrix failed!\n");
797
798 if (retries < 15)
799 {
800 if (retries == 0)
801 diagfact=0.1;
802 retries=retries+1;
803 diagfact=diagfact*10.0;
804 #pragma omp parallel for schedule(dynamic,64) private(i,j) shared(O,k,ldam)
805 for (i=1; i<=k-1; i++)
806 for (j=i; j<=k; j++)
807 O[ijtok(i,j, (long int) ldam)]=O[ijtok(j,i, (long int) ldam)];
808 for (i=1; i<=k; i++)
809 O[ijtok(i,i,(long int) ldam)]=diagO[i];
810 goto RETRYFACTOR;
811 }
812 else
813 {
814 if (printlevel >= 1)
815 printf("Factorization of the system matrix failed, giving up. \n");
816
817 /*
818 * Tighten up the solution as much as possible.
819 */
820
821 retcode=8;
822 goto RETURNBEST;
823
824 };
825 };
826
827
828 /*
829 Compute the rhs vector.
830 */
831
832 for (i=1; i<=k; i++)
833 rhs[i]=-a[i];
834
835 /*
836 * Add in the corrections for Fd.
837 */
838
839 /*
840 his section is where Fd is computed and put into
841 The rhs.
842
843 To save storage, work2 is used instead of a variable named
844 Fd.
845
846 dX will be used as a work variable where needed in place
847 of old work2 usage.
848 */
849
850 op_at(k,y,constraints,work1);
851
852 addscaledmat(Z,1.0,C,work2);
853
854 if ((bestmeas > 1.0e3) && (parameters.perturbobj>0))
855 {
856 if (printlevel >= 3)
857 printf("Perturbing C.\n");
858 make_i(work3);
859 addscaledmat(work2,-perturbfac,work3,work2);
860 };
861
862 if ((bestmeas < 1.0e3) && (parameters.perturbobj>0))
863 {
864 if (printlevel >= 3)
865 printf("Perturbing C.\n");
866 make_i(work3);
867 addscaledmat(work2,-perturbfac*pow(bestmeas/1000.0,1.5),work3,work2);
868 };
869
870
871
872 addscaledmat(work2,-1.0,work1,work2);
873
874 scale1=1.0;
875 scale2=0.0;
876
877
878
879 mat_multspb(scale1,scale2,Zi,work2,work3,fill);
880 mat_multspc(scale1,scale2,work3,X,dX,fill);
881
882
883
884 op_a(k,constraints,dX,workvec1);
885
886
887
888 for (i=1; i<=k; i++)
889 rhs[i] = rhs[i]+workvec1[i];
890
891 for (i=1; i<=k; i++)
892 workvec1[i]=rhs[i];
893
894 if (printlevel >=3)
895 {
896 printf("Fnorm of Fd is %e \n",Fnorm(work2));
897 printf("Norm of rhs is %e \n",norm2(m,rhs+1));
898 };
899
900 /*
901 Solve the system of equations for dy.
902 */
903
904 /*
905 * First, scale
906 */
907
908 for (i=1; i<=k; i++)
909 workvec1[i]=workvec1[i]*workvec8[i];
910
911 info=solvesys(k,ldam,O,workvec1);
912
913 for (i=1; i<=k; i++)
914 workvec1[i]=workvec1[i]*workvec8[i];
915
916 if (info != 0)
917 {
918 if (printlevel >= 1)
919 printf("Solving for dy failed! \n");
920 retcode=8;
921 goto RETURNBEST;
922 };
923
924 /*
925 * Do iterative refinement.
926 */
927
928
929 if ((iter>1) && (relerr2 > parameters.axtol) &&
930 (parameters.fastmode==0))
931 {
932 op_at(k,workvec1,constraints,work1);
933 mat_multspa(1.0,0.0,work1,X,dX,fill);
934 mat_multspc(1.0,0.0,Zi,dX,work1,fill);
935 op_a(k,constraints,work1,workvec2);
936 for (i=1; i<=k; i++)
937 workvec2[i]=rhs[i]-workvec2[i];
938
939 relerr=norm2(k,workvec2+1)/(1.0+norm2(k,rhs+1));
940 besterr=relerr;
941 for (i=1; i<=k; i++)
942 workvec4[i]=workvec1[i];
943
944 if (printlevel >= 3)
945 {
946 printf("refinement: Before relative error in Ody=r is %e \n",besterr);
947 fflush(stdout);
948 };
949
950 refinements=0;
951 lastimproverefinement=0;
952
953 while ((refinements < 20) && (refinements-lastimproverefinement < 3)
954 && (besterr > 1.0e-14))
955 {
956 refinements++;
957
958 for (i=1; i<=k; i++)
959 workvec3[i]=workvec2[i]*workvec8[i];
960
961 info=solvesys(k,ldam,O,workvec3);
962
963 for (i=1; i<=k; i++)
964 workvec3[i]=workvec3[i]*workvec8[i];
965
966 for (i=1; i<=k; i++)
967 workvec1[i]=workvec1[i]+workvec3[i];
968
969 op_at(k,workvec1,constraints,work1);
970 mat_multspa(1.0,0.0,work1,X,dX,fill);
971 mat_multspc(1.0,0.0,Zi,dX,work1,fill);
972 op_a(k,constraints,work1,workvec2);
973 for (i=1; i<=k; i++)
974 workvec2[i]=rhs[i]-workvec2[i];
975
976 relerr=norm2(k,workvec2+1)/(1.0+norm2(k,rhs+1));
977 if (relerr < besterr)
978 {
979 lastimproverefinement=refinements;
980 besterr=relerr;
981 for (i=1; i<=k; i++)
982 workvec4[i]=workvec1[i];
983 };
984 if (printlevel >= 4)
985 {
986 printf("refinement: During relative error in Ody=r is %e \n",besterr);
987 fflush(stdout);
988 };
989
990 };
991
992 if (printlevel >= 3)
993 printf("refinement: After relative error in Ody=r is %e, %d \n",besterr,lastimproverefinement);
994 for (i=1; i<=k; i++)
995 workvec1[i]=workvec4[i];
996
997 };
998
999
1000 /*
1001 * Extract dy.
1002 */
1003
1004 for (i=1; i<=k; i++)
1005 dy[i]=workvec1[i];
1006
1007 if (printlevel >= 3)
1008 printf("Norm of dy is %e \n",norm2(k,dy+1));
1009
1010 /*
1011 Compute dZ
1012 */
1013
1014 op_at(k,dy,constraints,dZ);
1015
1016 /*
1017 * Note: At this point, dZ only has A'(dy), not -Fd
1018 */
1019 /*
1020 Compute Zi*A'(dy), and save it in a temp variable for later use.
1021 */
1022
1023 scale1=1.0;
1024 scale2=0.0;
1025
1026 mat_multspb(scale1,scale2,Zi,dZ,work1,fill);
1027
1028 if (printlevel >= 4)
1029 printf("Fnorm of work1 is %e \n",Fnorm(work1));
1030
1031 /*
1032 * Now, update dZ to include -Fd
1033 */
1034
1035 if (printlevel >= 3)
1036 {
1037 printf("Before Fd, Fnorm of dZ is %e \n",Fnorm(dZ));
1038 fflush(stdout);
1039 };
1040
1041 addscaledmat(dZ,-1.0,work2,dZ);
1042
1043 if (printlevel >= 3)
1044 {
1045 printf("After Fd, Fnorm of dZ is %e \n",Fnorm(dZ));
1046 fflush(stdout);
1047 };
1048
1049
1050 /*
1051 * Now, we've got dZ in dZ, and Zi*A'(dy) in work1,
1052 * and Zi*Fd in work3
1053 */
1054
1055 /*
1056 Compute dX=-X+Zi*Fd*X-temp*X;
1057
1058 First, put I-Zi*Fd+work1 in workn2. Then multiply -work2*X, and
1059 put the result in dX.
1060
1061 */
1062 make_i(work2);
1063
1064 addscaledmat(work2,-1.0,work3,work2);
1065
1066 addscaledmat(work2,1.0,work1,work2);
1067
1068 scale1=-1.0;
1069 scale2=0.0;
1070 mat_mult(scale1,scale2,work2,X,dX);
1071
1072 sym_mat(dX);
1073
1074 if (printlevel >= 3)
1075 {
1076 printf("Fnorm of dX is %e \n",Fnorm(dX));
1077 printf("Fnorm of dZ is %e \n",Fnorm(dZ));
1078 fflush(stdout);
1079 };
1080
1081
1082
1083
1084 /*
1085 * Next, determine mu.
1086 */
1087
1088 if (relpinfeas < parameters.axtol)
1089 alphap1=linesearch(n,dX,work1,work2,work3,cholxinv,workvec4,
1090 workvec5,workvec6,mystepfrac,1.0,
1091 printlevel);
1092 else
1093 alphap1=linesearch(n,dX,work1,work2,work3,cholxinv,workvec4,
1094 workvec5,workvec6,mystepfrac,1.0,
1095 printlevel);
1096
1097 if (reldinfeas < parameters.atytol)
1098 alphad1=linesearch(n,dZ,work1,work2,work3,cholzinv,workvec4,
1099 workvec5,workvec6,mystepfrac,1.0,
1100 printlevel);
1101 else
1102 alphad1=linesearch(n,dZ,work1,work2,work3,cholzinv,workvec4,
1103 workvec5,workvec6,mystepfrac,1.0,
1104 printlevel);
1105
1106 /* Here, work1 holds X+alphap1*dX, work2=Z+alphad1*dZ */
1107
1108 addscaledmat(X,alphap1,dX,work1);
1109 addscaledmat(Z,alphad1,dZ,work2);
1110 for (i=1; i<=k; i++)
1111 workvec1[i]=y[i]+alphad1*dy[i];
1112
1113
1114 /*
1115 * Check to see whether this affine solution is the best yet.
1116 * The test is somewhat expensive, so don't do it unless
1117 * we're close to done.
1118 */
1119
1120 if ((bestmeas <1.0e4) && (parameters.affine==0))
1121 {
1122 /*
1123 * Verify that the new X and Z are Cholesky factorizable.
1124 */
1125
1126 copy_mat(work1,work3);
1127 ret=chol(work3);
1128 while (ret != 0)
1129 {
1130 if (printlevel >=3)
1131 printf("Affine eigsearch missed: adjusting alphap1\n");
1132 alphap1=alphap1*0.9;
1133 addscaledmat(X,alphap1,dX,work1);
1134 copy_mat(work1,work3);
1135 ret=chol(work3);
1136 };
1137
1138 copy_mat(work2,work3);
1139 ret=chol(work3);
1140 while (ret != 0)
1141 {
1142 if (printlevel >=3)
1143 printf("Affine eigsearch missed: adjusting alphad1\n");
1144 alphad1=alphad1*0.9;
1145 addscaledmat(Z,alphad1,dZ,work2);
1146 for (i=1; i<=k; i++)
1147 workvec1[i]=y[i]+alphad1*dy[i];
1148 copy_mat(work2,work3);
1149 ret=chol(work3);
1150 };
1151
1152 /*
1153 * Now, check the quality of this solution.
1154 */
1155
1156 affpobj=calc_pobj(C,work1,constant_offset);
1157 affdobj=calc_dobj(k,a,workvec1,constant_offset);
1158
1159 /*
1160 * run user exit to check if the affine solution is good enough
1161 */
1162 ret=user_exit(n,k,C,a,affdobj,affpobj,constant_offset,constraints,
1163 work1,workvec1,work2,parameters);
1164 if (ret >= 1)
1165 {
1166 if (ret == 1)
1167 {
1168 /*
1169 * Received SIGTERM, SIGQUIT, ...
1170 */
1171 return(10);
1172 }
1173 else
1174 {
1175 *dobj=affdobj;
1176 *pobj=affpobj;
1177 copy_mat(work1,X);
1178 copy_mat(work2,Z);
1179 for(i=1; i<=k; i++)
1180 y[i]=workvec1[i];
1181 if(printlevel>=3)
1182 printf("Affine step good enough, exiting\n");
1183 return(0);
1184 };
1185 };
1186
1187 if (parameters.usexzgap==0)
1188 {
1189 affgap=affdobj-affpobj;
1190 if (affgap < 0)
1191 affgap=0.0;
1192 affrelgap=affgap/(1.0+fabs(affpobj)+fabs(affdobj));
1193 }
1194 else
1195 {
1196 affgap=trace_prod(work1,work2);
1197 if (affgap < 0)
1198 affgap=0.0;
1199 affrelgap=affgap/(1.0+fabs(affpobj)+fabs(affdobj));
1200 };
1201
1202 affreldinfeas=dinfeas(k,C,constraints,workvec1,work2,work3);
1203 affrelpinfeas=pinfeas(k,constraints,work1,a,workvec4);
1204
1205 if (printlevel >= 3)
1206 {
1207 printf("affpobj is %e \n",affpobj);
1208 printf("affdobj is %e \n",affdobj);
1209 printf("affrelgap is %e \n",affrelgap);
1210 printf("affrelpinfeas is %e \n",affrelpinfeas);
1211 printf("affreldinfeas is %e \n",affreldinfeas);
1212 };
1213
1214 if ((affrelgap/parameters.objtol < bestmeas) &&
1215 (affrelpinfeas/parameters.axtol <bestmeas) &&
1216 (affreldinfeas/parameters.atytol <bestmeas))
1217 {
1218 bestmeas=affrelgap/parameters.objtol;
1219 if (affrelpinfeas/parameters.axtol > bestmeas)
1220 bestmeas=affrelpinfeas/parameters.axtol;
1221 if (affreldinfeas/parameters.atytol > bestmeas)
1222 bestmeas=affreldinfeas/parameters.atytol;
1223
1224 store_packed(work1,bestx);
1225 store_packed(work2,bestz);
1226 for (i=1; i<=k; i++)
1227 besty[i]=y[i]+alphad1*dy[i];
1228 if (printlevel >= 3)
1229 printf("Affine step: New best solution, %e \n",bestmeas);
1230 };
1231
1232 if ((ispfeasprob==1) &&
1233 (affrelpinfeas/parameters.axtol < bestmeas))
1234 {
1235 bestmeas=affrelpinfeas/parameters.axtol;
1236
1237 store_packed(work1,bestx);
1238 for (i=1; i<=k; i++)
1239 besty[i]=0.0;
1240 zero_mat(work3);
1241 addscaledmat(work3,1.0e-50,work1,work3);
1242 store_packed(work3,bestz);
1243
1244 if (printlevel >= 3)
1245 printf("Affine step: New best solution, %e \n",bestmeas);
1246 };
1247
1248 if ((isdfeasprob==1) &&
1249 (affreldinfeas/parameters.atytol <bestmeas))
1250 {
1251 bestmeas=affreldinfeas/parameters.atytol;
1252
1253 zero_mat(work3);
1254 make_i(work1);
1255 addscaledmat(work3,1.0e-40,work1,work3);
1256 store_packed(work3,bestx);
1257 store_packed(work2,bestz);
1258 for (i=1; i<=k; i++)
1259 besty[i]=workvec1[i];
1260 if (printlevel >= 3)
1261 printf("Affine step: New best solution, %e \n",bestmeas);
1262 };
1263
1264 if (bestmeas < 1.0)
1265 {
1266 if (printlevel >= 3)
1267 printf("Finishing with a final affine step.\n");
1268 iter=iter+1;
1269 if (printlevel >= 1)
1270 printf("Iter: %2d Ap: %.2e Pobj: % .7e Ad: %.2e Dobj: % .7e \n",iter,alphap1,affpobj,alphad1,affdobj);
1271 if (printlevel >= 2)
1272 printf("Total Iterations: %d \n",iter);
1273
1274 store_unpacked(bestx,X);
1275 store_unpacked(bestz,Z);
1276 for (i=1; i<=k; i++)
1277 y[i]=besty[i];
1278 *pobj=calc_pobj(C,X,constant_offset);
1279 *dobj=calc_dobj(k,a,y,constant_offset);
1280 return(0);
1281 };
1282 };
1283
1284 /*
1285 * Compute muplus and prepare for the corrector step.
1286 */
1287
1288 muplus=trace_prod(work1,work2)/(n);
1289 muk=trace_prod(X,Z)/(n);
1290
1291 if (muk < 0.0)
1292 muk=fabs(muk);
1293 if (muplus < 0.0)
1294 muplus=muk/2;
1295
1296 gamma=(muplus/muk);
1297
1298 /*
1299 * Pick the new mu as follows:
1300 *
1301 * If we have been making good progress (alphap > 0.5) and
1302 * (alphad>0.5) then use mu=muk*gamma*gamma*min(gamma,0.5);
1303 *
1304 * If we haven't been making good progress, then just do
1305 * mu=muplus.
1306 *
1307 * Also, make sure that mu is no larger than the old mu.
1308 */
1309
1310 if ((relpinfeas < 0.1*parameters.axtol) && (alphad>0.2) &&
1311 (reldinfeas < 0.1*parameters.atytol) && (alphap>0.2) &&
1312 (mu > 1.0e-6) && (gamma < 1.0) && (alphap+alphad > 1.0))
1313 {
1314 mu=muk*pow(gamma,alphap+alphad);
1315 }
1316 else
1317 {
1318 if (muplus < 0.9*muk)
1319 mu=muplus;
1320 else
1321 mu=muk*0.9;
1322 };
1323
1324 /*
1325 * If we have primal and dual infeasibility in hand, then
1326 * make sure that mu is <=muk/2.
1327 */
1328
1329 if ((relpinfeas < 0.9*parameters.axtol) &&
1330 (reldinfeas < 0.9*parameters.atytol) &&
1331 (mu > muk/2))
1332 mu=muk/2;
1333
1334 /*
1335 * If we want a primal-dual affine step, then set mu=0.0.
1336 */
1337
1338 if (parameters.affine==1)
1339 {
1340 mu=0.0;
1341 if (printlevel >= 3)
1342 printf("Taking an affine step because parameters.affine=1\n");
1343 };
1344
1345
1346 /*
1347 * Printout some info on mu.
1348 */
1349
1350 if (printlevel >= 2)
1351 {
1352 printf("muk is %e \n",muk);
1353 printf("muplus is %e \n",muplus);
1354 printf("New mu is %e \n",mu);
1355 printf("mu*n = target duality gap is %e \n",mu*n);
1356 fflush(stdout);
1357 };
1358
1359 /*
1360 * Take a moment to figure out how well we're doing on feasibility.
1361 */
1362
1363 addscaledmat(X,1.0,dX,work1);
1364 op_a(k,constraints,work1,workvec1);
1365 for (i=1; i<=k; i++)
1366 workvec1[i]=workvec1[i]-a[i];
1367 relerr1=norm2(k,workvec1+1)/(1.0+norma);
1368 if (printlevel >= 3)
1369 {
1370 printf("refinement: Relative error in A(X+dX)=a (Fphat) is %e \n",
1371 relerr1);
1372 };
1373
1374 /*
1375 Now, compute the corrector step.
1376 */
1377
1378 /*
1379 rhs for the corrector step.
1380 */
1381
1382 /*
1383 Update Fphat=a-A(X+dX)
1384 */
1385
1386 addscaledmat(dX,1.0,X,work1);
1387 op_a(k,constraints,work1,Fp);
1388 for (i=1; i<=k; i++)
1389 Fp[i]=a[i]-Fp[i];
1390
1391 /*
1392 The RHS is now A(Zi*Fdhat*X)+A(Zi*(-dZ*dX+mu*I))-Fphat
1393
1394 = A(Zi*(Fdhat*X-dZ*dX+mu*I))-Fphat
1395
1396 */
1397
1398 make_i(work1);
1399 scale1=0.0;
1400 scale2=mu;
1401 mat_mult(scale1,scale2,work2,work2,work1);
1402
1403 scale1=-1.0;
1404 scale2=1.0;
1405 mat_multspa(scale1,scale2,dZ,dX,work1,fill);
1406
1407 scale1=1.0;
1408 scale2=0.0;
1409 mat_multspc(scale1,scale2,Zi,work1,work2,fill);
1410
1411
1412
1413 /*
1414 Next, compute op_a of work2, and put the result in rhs.
1415 */
1416
1417 op_a(k,constraints,work2,rhs);
1418
1419 /*
1420 * Finally, subtract off Fphat.
1421 */
1422
1423 for (i=1; i<=k; i++)
1424 rhs[i]=rhs[i]-Fp[i];
1425
1426 for (i=1; i<=k; i++)
1427 workvec1[i]=rhs[i];
1428
1429 /*
1430 Solve for dy1.
1431 */
1432
1433 for (i=1; i<=k; i++)
1434 workvec1[i]=workvec1[i]*workvec8[i];
1435
1436 info=solvesys(k,ldam,O,workvec1);
1437
1438 for (i=1; i<=k; i++)
1439 workvec1[i]=workvec1[i]*workvec8[i];
1440
1441 if (info != 0)
1442 {
1443 if (printlevel >= 1)
1444 printf("Solving for dy1 failed! \n");
1445 retcode=8;
1446 goto RETURNBEST;
1447 };
1448
1449 /*
1450 * Do iterative refinement.
1451 */
1452
1453 if ((iter>1) && (relerr2 > 0.01*parameters.axtol) &&
1454 (parameters.fastmode==0))
1455 {
1456 op_at(k,workvec1,constraints,work1);
1457 mat_multspa(1.0,0.0,work1,X,work2,fill);
1458 mat_multspc(1.0,0.0,Zi,work2,work3,fill);
1459 op_a(k,constraints,work3,workvec2);
1460 for (i=1; i<=k; i++)
1461 workvec2[i]=rhs[i]-workvec2[i];
1462
1463 relerr=norm2(k,workvec2+1)/(1.0+norm2(k,rhs+1));
1464 besterr=relerr;
1465 for (i=1; i<=k; i++)
1466 workvec4[i]=workvec1[i];
1467
1468 if (printlevel >= 3)
1469 {
1470 printf("refinement: Before relative error in Odyhat=r is %e \n",besterr);
1471 };
1472
1473 refinements=0;
1474 lastimproverefinement=0;
1475
1476 while ((refinements < 20) && (refinements-lastimproverefinement < 3)
1477 && (besterr > 1.0e-14))
1478 {
1479 refinements++;
1480
1481 for (i=1; i<=k; i++)
1482 workvec3[i]=workvec2[i]*workvec8[i];
1483
1484 info=solvesys(k,ldam,O,workvec3);
1485
1486 for (i=1; i<=k; i++)
1487 workvec3[i]=workvec3[i]*workvec8[i];
1488
1489 for (i=1; i<=k; i++)
1490 workvec1[i]=workvec1[i]+workvec3[i];
1491
1492 op_at(k,workvec1,constraints,work1);
1493 mat_multspa(1.0,0.0,work1,X,work2,fill);
1494 mat_multspc(1.0,0.0,Zi,work2,work3,fill);
1495 op_a(k,constraints,work3,workvec2);
1496 for (i=1; i<=k; i++)
1497 workvec2[i]=rhs[i]-workvec2[i];
1498
1499 relerr=norm2(k,workvec2+1)/(1.0+norm2(k,rhs+1));
1500 if (relerr < besterr)
1501 {
1502 lastimproverefinement=refinements;
1503 besterr=relerr;
1504 for (i=1; i<=k; i++)
1505 workvec4[i]=workvec1[i];
1506 };
1507 if (printlevel >= 4)
1508 {
1509 printf("refinement: During relative error in Ody=r is %e \n",besterr);
1510 fflush(stdout);
1511 };
1512
1513 };
1514
1515 if (printlevel >= 3)
1516 {
1517 printf("refinement: After relative error in Odyhat=r is %e,%d \n",besterr,lastimproverefinement);
1518 };
1519
1520 for (i=1; i<=k; i++)
1521 workvec1[i]=workvec4[i];
1522 };
1523
1524 /*
1525 * retrieve dy1.
1526 */
1527
1528 for (i=1; i<=k; i++)
1529 dy1[i]=workvec1[i];
1530
1531 /*
1532 Compute dZ1=A'(dy1).
1533
1534 dZ1 is stored in work3.
1535 */
1536
1537 op_at(k,dy1,constraints,work3);
1538
1539
1540 /*
1541 Compute dX1=-Zi*dZ1*X-Zi*dZ*dX+mu*zi;
1542 dX1=Zi*(-dZ1*X-dZ*dX+mu*I)
1543
1544 for storage sake, dX1 is stored in work2.
1545 */
1546
1547 make_i(work1);
1548 scale1=-1.0;
1549 scale2=mu;
1550 mat_multspa(scale1,scale2,dZ,dX,work1,fill);
1551 scale1=-1.0;
1552 scale2=1.0;
1553 mat_multspa(scale1,scale2,work3,X,work1,fill);
1554 scale1=1.0;
1555 scale2=0.0;
1556 mat_mult(scale1,scale2,Zi,work1,work2);
1557 sym_mat(work2);
1558
1559 addscaledmat(X,1.0,dX,work1);
1560 op_a(k,constraints,work1,workvec1);
1561 for (i=1; i<=k; i++)
1562 workvec1[i]=workvec1[i]-a[i];
1563 relerr2=norm2(k,workvec1+1)/(1.0+norma);
1564
1565 if (printlevel >= 3)
1566 {
1567 printf("refinement: Before dX+dX1 Relative error in A(X+dX)=a is %e \n",relerr2);
1568 if (relerr1 < relerr2)
1569 printf("refinement: worse\n");
1570 else
1571 printf("refinement: better\n");
1572 };
1573 /*
1574 Update the predictor step.
1575 */
1576
1577 if (printlevel >= 3)
1578 {
1579 printf("Fnorm of dX1 is %e\n",Fnorm(work2));
1580 printf("Fnorm of dZ1 is %e\n",Fnorm(work3));
1581 };
1582 add_mat(work2,dX);
1583 add_mat(work3,dZ);
1584 for (i=1; i<=k; i++)
1585 dy[i]=dy1[i]+dy[i];
1586
1587
1588
1589 /*
1590 * Check A(X+dX)=a.
1591 */
1592 addscaledmat(X,1.0,dX,work1);
1593 op_a(k,constraints,work1,workvec1);
1594 for (i=1; i<=k; i++)
1595 workvec1[i]=workvec1[i]-a[i];
1596 relerr2=norm2(k,workvec1+1)/(1.0+norma);
1597
1598 if (printlevel >= 3)
1599 {
1600 printf("refinement: Before adjust dX Relative error in A(X+dX)=a is %e \n",relerr2);
1601 if (relerr1 < relerr2)
1602 printf("refinement: worse\n");
1603 else
1604 printf("refinement: better\n");
1605 };
1606
1607 if (printlevel >= 3)
1608 {
1609 printf("Fnorm of dX is %e \n",Fnorm(dX));
1610 printf("Fnorm of dZ is %e \n",Fnorm(dZ));
1611 };
1612
1613 /*
1614 * Take a moment to figure out how well we're doing on feasibility.
1615 */
1616
1617 addscaledmat(X,1.0,dX,work1);
1618 op_a(k,constraints,work1,workvec1);
1619 for (i=1; i<=k; i++)
1620 workvec1[i]=workvec1[i]-a[i];
1621 relerr2=norm2(k,workvec1+1)/(1.0+norma);
1622
1623 if (printlevel >= 3)
1624 {
1625 printf("refinement: After adjust error in A(X+dX)=a is %e \n",relerr2);
1626 if (relerr1 < relerr2)
1627 printf("refinement: worse\n");
1628 else
1629 printf("refinement: better\n");
1630 };
1631
1632
1633 /*
1634 Now, we've got the individual steps.
1635 Find maximum possible step sizes.
1636 */
1637
1638 alphap=linesearch(n,dX,work1,work2,work3,cholxinv,workvec4,
1639 workvec5,workvec6,mystepfrac,1.0,
1640 printlevel);
1641
1642 alphad=linesearch(n,dZ,work1,work2,work3,cholzinv,workvec4,
1643 workvec5,workvec6,mystepfrac,1.0,
1644 printlevel);
1645
1646 if (printlevel >= 3)
1647 {
1648 printf("After linesearches, alphap=%e alphad=%e \n",alphap,alphad);
1649 };
1650
1651 for (i=1; i<=k; i++)
1652 workvec1[i]=y[i]+alphad*dy[i];
1653
1654 /*
1655 * Check on the feasibility of the new solutions. If they're
1656 * worse, and the old solution was feasible, then don't take the
1657 * step.
1658 */
1659
1660 addscaledmat(X,alphap,dX,work1);
1661 newrelpinfeas=pinfeas(k,constraints,work1,a,workvec1);
1662
1663 /*
1664 * For the primal infeasibility, check the relative gap and the
1665 * current primal infeasibility to establish a limit on the
1666 * new infeasibility.
1667 */
1668
1669 limrelpinfeas=bestrelpinfeas*100;
1670
1671 if ((limrelpinfeas >relgap) && (relgap > 0))
1672 limrelpinfeas=relgap;
1673
1674 /*
1675 * In the early stages, don't worry about limrelpinfeas.
1676 */
1677
1678 if ((iter < 10) || (relpinfeas > 0.01) ||
1679 (dinfeasmeas > 1.0e-3*parameters.dinftol))
1680 limrelpinfeas=relpinfeas*1000;
1681
1682 /*
1683 * Don't ever ask for more than the required tolerance.
1684 */
1685
1686 if (parameters.axtol > limrelpinfeas)
1687 limrelpinfeas=parameters.axtol;
1688
1689 /*
1690 * If we're trying to prove dual infeasibility than don't limit
1691 * the step.
1692 */
1693
1694 if ((probdfeas==0) && (dinfeasmeas>1.0e4))
1695 limrelpinfeas=1.0e30;
1696
1697 /*
1698 * Now, make sure that the step keeps us feasible enough.
1699 */
1700
1701 if (printlevel >= 3)
1702 {
1703 printf("newrelpinfeas is %e\n",newrelpinfeas);
1704 printf("limrelpinfeas is %e\n",limrelpinfeas);
1705 };
1706
1707 i=1;
1708 while (newrelpinfeas > limrelpinfeas)
1709 {
1710 alphap=0.80*alphap;
1711 addscaledmat(X,alphap,dX,work1);
1712 newrelpinfeas=pinfeas(k,constraints,work1,a,workvec1);
1713 i=i+1;
1714
1715 if (i>20)
1716 {
1717 if (printlevel >=3)
1718 printf("Stuck at edge of primal feasibility.\n");
1719
1720 if (retries < 15)
1721 {
1722 if (retries == 0)
1723 diagfact=0.1;
1724 retries=retries+1;
1725 diagfact=diagfact*10.0;
1726 for (i=1; i<=k-1; i++)
1727 for (j=i; j<=k; j++)
1728 O[ijtok(i,j, (long int) ldam)]=O[ijtok(j,i, (long int) ldam)];
1729 for (i=1; i<=k; i++)
1730 O[ijtok(i,i,(long int) ldam)]=diagO[i];
1731 goto RETRYFACTOR;
1732 }
1733 else
1734 {
1735 if (printlevel >= 1)
1736 printf("Stuck at edge of primal feasibility, giving up. \n");
1737
1738 /*
1739 * Tighten up the solution as much as possible.
1740 */
1741
1742 retcode=5;
1743 goto RETURNBEST;
1744
1745 };
1746
1747
1748 };
1749 };
1750
1751 /*
1752 Now make sure that we aren't stepping outside of dual feasibility.
1753 */
1754
1755
1756 addscaledmat(Z,alphad,dZ,work1);
1757
1758 for (i=1; i<=k; i++)
1759 workvec1[i]=y[i]+alphad*dy[i];
1760
1761 newreldinfeas=dinfeas(k,C,constraints,workvec1,work1,work2);
1762
1763 limreldinfeas=bestreldinfeas*100;
1764
1765 if ((limreldinfeas > relgap) && (relgap > 0))
1766 limreldinfeas=relgap;
1767
1768 /*
1769 * In the early stages, don't worry about limreldinfeas.
1770 * Also don't worry if we're trying to establish primal
1771 * infeasibility.
1772 */
1773
1774 if ((iter < 10) || (reldinfeas > 0.01) ||
1775 (pinfeasmeas > 1.0e-3*parameters.pinftol))
1776 limreldinfeas=reldinfeas*1000;
1777
1778 /*
1779 * Don't ever ask for more than the required tolerance.
1780 */
1781
1782 if (parameters.atytol > limreldinfeas)
1783 limreldinfeas=parameters.atytol;
1784
1785 if (printlevel >= 3)
1786 {
1787 printf("newreldinfeas is %e\n",newreldinfeas);
1788 printf("limreldinfeas is %e\n",limreldinfeas);
1789 };
1790
1791 i=1;
1792 while (newreldinfeas > limreldinfeas)
1793 {
1794 alphad=0.80*alphad;
1795 addscaledmat(Z,alphad,dZ,work1);
1796
1797 for (j=1; j<=k; j++)
1798 workvec1[j]=y[j]+alphad*dy[j];
1799
1800 newreldinfeas=dinfeas(k,C,constraints,workvec1,work1,
1801 work2);
1802 i=i+1;
1803 if (i>15)
1804 {
1805 if (printlevel >=3)
1806 printf("Stuck at edge of dual feasibility.\n");
1807
1808 if (retries < 15)
1809 {
1810 if (retries == 0)
1811 diagfact=0.1;
1812 retries=retries+1;
1813 diagfact=diagfact*10.0;
1814 for (i=1; i<=k-1; i++)
1815 for (j=i; j<=k; j++)
1816 O[ijtok(i,j,(long int) ldam)]=O[ijtok(j,i,(long int) ldam)];
1817 for (i=1; i<=k; i++)
1818 O[ijtok(i,i,(long int) ldam)]=diagO[i];
1819 goto RETRYFACTOR;
1820 }
1821 else
1822 {
1823 if (printlevel >= 1)
1824 printf("Stuck at edge of dual feasibility, giving up. \n");
1825
1826 /*
1827 * Tighten up the solution as much as possible.
1828 */
1829
1830 retcode=6;
1831 goto RETURNBEST;
1832
1833 };
1834
1835 };
1836 };
1837
1838 if (printlevel >= 3)
1839 {
1840 printf("After feas check, alphap=%e alphad=%e \n",alphap,alphad);
1841 };
1842
1843 /*
1844 * Give up if step lengths are way too small.
1845 */
1846
1847 if ((alphap<=parameters.minstepp) ||
1848 (alphad<=parameters.minstepd))
1849 {
1850 if (printlevel >= 2)
1851 printf("line search failure in corrector step.\n");
1852
1853 if (retries < 15)
1854 {
1855 if (retries == 0)
1856 diagfact=0.1;
1857 retries=retries+1;
1858 diagfact=diagfact*10.0;
1859 for (i=1; i<=k-1; i++)
1860 for (j=i; j<=k; j++)
1861 O[ijtok(i,j,(long int) ldam)]=O[ijtok(j,i,(long int) ldam)];
1862 for (i=1; i<=k; i++)
1863 O[ijtok(i,i,(long int) ldam)]=diagO[i];
1864 goto RETRYFACTOR;
1865 }
1866 else
1867 {
1868 if (printlevel >= 1)
1869 printf("Too many line search failures, giving up. \n");
1870
1871 /*
1872 * Tighten up the solution as much as possible.
1873 */
1874
1875 retcode=7;
1876 goto RETURNBEST;
1877
1878 };
1879 };
1880
1881 /*
1882 * In case alphap changed, recompute these.
1883 */
1884
1885 addscaledmat(X,alphap,dX,work1);
1886 newrelpinfeas=pinfeas(k,constraints,work1,a,workvec1);
1887
1888 addscaledmat(Z,alphad,dZ,work1);
1889 for (i=1; i<=k; i++)
1890 workvec1[i]=y[i]+alphad*dy[i];
1891
1892 newreldinfeas=dinfeas(k,C,constraints,workvec1,work1,work2);
1893
1894 /*
1895 * Update cholzinv.
1896 */
1897
1898 copy_mat(work1,work2);
1899 ret=chol(work2);
1900
1901 while (ret != 0)
1902 {
1903 alphad=alphad*0.90;
1904 addscaledmat(Z,alphad,dZ,work1);
1905 for (i=1; i<=k; i++)
1906 workvec1[i]=y[i]+alphad*dy[i];
1907 copy_mat(work1,work2);
1908 ret=chol(work2);
1909 if (printlevel >=2)
1910 printf("eigsearch missed! Adjusting alphad\n");
1911 };
1912
1913 chol_inv(work2,work3);
1914 store_packed(work3,cholzinv);
1915
1916 /*
1917 Compute Zi.
1918 */
1919
1920 copy_mat(work3,work1);
1921 trans(work1);
1922 scale1=1.0;
1923 scale2=0.0;
1924 mat_mult(scale1,scale2,work3,work1,Zi);
1925 if (printlevel >= 4)
1926 printf("Fnorm of Zi is %e \n",Fnorm(Zi));
1927
1928 /*
1929 * Confirm that X is in fact factorable. If not, reduce
1930 * alpha until it is.
1931 */
1932
1933 addscaledmat(X,alphap,dX,work1);
1934 copy_mat(work1,work2);
1935 ret=chol(work2);
1936
1937 while (ret != 0)
1938 {
1939 if (printlevel >=2)
1940 printf("eigsearch missed! Adjusting alphap\n");
1941 alphap=alphap*0.90;
1942 addscaledmat(X,alphap,dX,work1);
1943 copy_mat(work1,work2);
1944 ret=chol(work2);
1945 };
1946
1947 chol_inv(work2,work3);
1948 store_packed(work3,cholxinv);
1949
1950 /*
1951 * do a line search for feasibility.
1952 */
1953
1954 if (relpinfeas > 1.0)
1955 {
1956 newalphap=alphap;
1957 for (i=1; i<=9; i++)
1958 {
1959 addscaledmat(X,(i*1.0)*alphap/10.0,dX,work1);
1960 if (pinfeas(k,constraints,work1,a,workvec1) < newrelpinfeas)
1961 {
1962 newalphap=i*1.0*alphap/10.0;
1963 newrelpinfeas=pinfeas(k,constraints,work1,a,workvec1);
1964 };
1965 };
1966 if (newalphap < alphap)
1967 {
1968 if (printlevel >= 2)
1969 printf("Feasibility Adjusting alphap to %e \n",newalphap);
1970 alphap=newalphap;
1971 };
1972 };
1973
1974 /*
1975 *Take the step.
1976 */
1977
1978 addscaledmat(X,alphap,dX,X);
1979 addscaledmat(Z,alphad,dZ,Z);
1980
1981 for (i=1; i<=k; i++)
1982 y[i]=y[i]+alphad*dy[i];
1983
1984 /*
1985 Recompute the objective function values.
1986 */
1987
1988 *pobj=calc_pobj(C,X,constant_offset);
1989 *dobj=calc_dobj(k,a,y,constant_offset);
1990 if (parameters.usexzgap==0)
1991 {
1992 gap=*dobj-*pobj;
1993 if (gap < 0.0)
1994 gap=0.0;
1995 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
1996 }
1997 else
1998 {
1999 gap=trace_prod(X,Z);
2000 if (gap < 0.0)
2001 gap=0.0;
2002 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
2003 };
2004
2005 if (printlevel >= 2)
2006 {
2007 printf("pobj is %e \n",*pobj);
2008 printf("dobj is %e \n",*dobj);
2009 printf("gap is %e \n",gap);
2010 printf("relgap is %e \n",relgap);
2011 };
2012
2013 relpinfeas=pinfeas(k,constraints,X,a,workvec1);
2014 reldinfeas=dinfeas(k,C,constraints,y,Z,work2);
2015
2016 if (relpinfeas < parameters.axtol )
2017 probpfeas=1;
2018
2019 if (reldinfeas < parameters.atytol)
2020 probdfeas=1;
2021
2022 /*
2023 * Make sure that the objective value hasn't gone crazy.
2024 *
2025 * This was a test with isnan, but isnan isn't ANSI C, so
2026 * we use an equivalent that typically works.
2027 *
2028 */
2029
2030 if ((gap) != (gap))
2031 {
2032 retcode=9;
2033 goto RETURNBEST;
2034 };
2035
2036 /*
2037 * This was a test with isinf, but isinf isn't ANSI C, so
2038 * we just test for extremely large values.
2039 */
2040
2041 if ((gap > 1.0e100) || (gap < -1.0e100))
2042 {
2043 retcode=9;
2044 goto RETURNBEST;
2045 };
2046
2047 /*
2048 * if this solution is better than the previous best, then
2049 * update our best solution.
2050 */
2051
2052 if ((relgap/parameters.objtol < bestmeas) &&
2053 (relpinfeas/parameters.axtol <bestmeas) &&
2054 (ispfeasprob==0) && (isdfeasprob==0) &&
2055 (reldinfeas/parameters.atytol <bestmeas))
2056 {
2057 bestmeas=relgap/parameters.objtol;
2058 if (relpinfeas/parameters.axtol > bestmeas)
2059 bestmeas=relpinfeas/parameters.axtol;
2060 if (reldinfeas/parameters.atytol > bestmeas)
2061 bestmeas=reldinfeas/parameters.atytol;
2062
2063 store_packed(X,bestx);
2064 store_packed(Z,bestz);
2065 for (i=1; i<=k; i++)
2066 besty[i]=y[i];
2067 if (printlevel >= 3)
2068 printf("New best solution, %e \n",bestmeas);
2069 };
2070
2071 if ((ispfeasprob==1) &&
2072 (relpinfeas/parameters.axtol <bestmeas))
2073 {
2074 bestmeas=relpinfeas/parameters.axtol;
2075
2076 store_packed(X,bestx);
2077 for (i=1; i<=k; i++)
2078 besty[i]=0.0;
2079 make_i(work2);
2080 zero_mat(work3);
2081 addscaledmat(work3,0.005*parameters.objtol/trace_prod(X,work2),work2,work3);
2082 store_packed(work3,bestz);
2083 if (printlevel >= 3)
2084 printf("New best solution, %e \n",bestmeas);
2085 };
2086
2087 if ((isdfeasprob==1) &&
2088 (reldinfeas/parameters.atytol <bestmeas))
2089 {
2090 bestmeas=reldinfeas/parameters.atytol;
2091
2092 zero_mat(work3);
2093 make_i(work1);
2094 addscaledmat(work3,1.0e-40,work1,work3);
2095 store_packed(work3,bestx);
2096 store_packed(Z,bestz);
2097 for (i=1; i<=k; i++)
2098 besty[i]=y[i];
2099 if (printlevel >= 3)
2100 printf("New best solution, %e \n",bestmeas);
2101 };
2102
2103 /*
2104 * If we haven't improved the gap much in 20 iterations, then
2105 * give up, because we're not making progress.
2106 */
2107
2108 if ((iter > 60) && (bestmeas > 0.5*bestarray[((iter-20) % BASIZE)]))
2109 {
2110 if (printlevel >= 1)
2111 printf("Lack of progress. Giving up!\n");
2112 retcode=7;
2113 goto RETURNBEST;
2114 };
2115
2116 /*
2117 Update counters and display status.
2118 */
2119
2120 iter++;
2121 if (printlevel >= 1)
2122 {
2123 printf("Iter: %2d Ap: %.2e Pobj: % .7e Ad: %.2e Dobj: % .7e \n",iter,alphap,*pobj,alphad,*dobj);
2124 fflush(stdout);
2125 };
2126
2127 /*
2128 * If iter gets above maxiter, then exit.
2129 */
2130
2131 if (iter >= parameters.maxiter)
2132 {
2133 if (printlevel >= 1)
2134 printf("Maximum iterations reached. \n");
2135 retcode=4;
2136 goto RETURNBEST;
2137 };
2138
2139 };
2140
2141 /*
2142 * Return success.
2143 */
2144
2145 retcode=0;
2146
2147 RETURNBEST:
2148
2149 store_unpacked(bestx,X);
2150 store_unpacked(bestz,Z);
2151 for (i=1; i<=k; i++)
2152 y[i]=besty[i];
2153
2154 if (parameters.usexzgap==0)
2155 {
2156 *pobj=calc_pobj(C,X,constant_offset);
2157 *dobj=calc_dobj(k,a,y,constant_offset);
2158 gap=*dobj-*pobj;
2159 }
2160 else
2161 {
2162 *pobj=calc_pobj(C,X,constant_offset);
2163 *dobj=calc_dobj(k,a,y,constant_offset);
2164 gap=trace_prod(X,Z);
2165 };
2166
2167
2168 if ((gap < 0.0) && (parameters.usexzgap==0) && (parameters.tweakgap==1)
2169 && (ispfeasprob==0) && (isdfeasprob==0))
2170 {
2171 tweakgap(n,k,a,constraints,gap,Z,dZ,y,dy,work1,work2,work3,dX,
2172 workvec1,workvec2,workvec3,workvec4,printlevel);
2173 };
2174
2175 /*
2176 Recompute the objective function values.
2177 */
2178
2179 *pobj=calc_pobj(C,X,constant_offset);
2180 *dobj=calc_dobj(k,a,y,constant_offset);
2181
2182 if (parameters.usexzgap==0)
2183 {
2184 gap=*dobj-*pobj;
2185 if (gap < 0.0)
2186 gap=0.0;
2187 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
2188 }
2189 else
2190 {
2191 gap=trace_prod(X,Z);
2192 if (gap < 0.0)
2193 gap=0.0;
2194 relgap=gap/(1.0+fabs(*dobj)+fabs(*pobj));
2195 };
2196
2197 /*
2198 * Recheck the primal and dual infeasibilities and gap.
2199 */
2200
2201 relpinfeas=pinfeas(k,constraints,X,a,workvec1);
2202 reldinfeas=dinfeas(k,C,constraints,y,Z,work1);
2203
2204 if (relpinfeas < parameters.axtol )
2205 probpfeas=1;
2206
2207 if (reldinfeas < parameters.atytol)
2208 probdfeas=1;
2209
2210 if ((relgap < parameters.objtol) &&
2211 (relpinfeas <parameters.axtol ) &&
2212 (reldinfeas <parameters.atytol))
2213 {
2214 /*
2215 * Our best solution actually did satisfy the conditions, so we've
2216 * succeeded in solving the problem!
2217 *
2218 */
2219 retcode=0;
2220 };
2221
2222 /*
2223 * Check for a good solution with slightly reduced accuracy.
2224 */
2225
2226 if ((ispfeasprob==0) && (isdfeasprob==0))
2227 {
2228 bestmeas=relgap/parameters.objtol;
2229 if ((relpinfeas/parameters.axtol) > bestmeas)
2230 bestmeas=relpinfeas/parameters.axtol;
2231 if ((reldinfeas/parameters.atytol) > bestmeas)
2232 bestmeas=reldinfeas/parameters.atytol;
2233
2234 if ((bestmeas > 1.0) && (bestmeas < 1000.0))
2235 retcode=3;
2236 };
2237
2238 if (ispfeasprob==1)
2239 {
2240 bestmeas=relpinfeas/parameters.axtol;
2241
2242 if ((bestmeas > 1.0) && (bestmeas < 1000.0))
2243 retcode=3;
2244 };
2245
2246 if (isdfeasprob==1)
2247 {
2248 bestmeas=reldinfeas/parameters.atytol;
2249
2250 if ((bestmeas > 1.0) && (bestmeas < 1000.0))
2251 retcode=3;
2252 };
2253
2254 /*
2255 * Come here if we have an infeasible problem and are returning the
2256 * certificate of infeasibility.
2257 */
2258
2259 RETURNCERT:
2260
2261 /*
2262 * Print out the total number of iterations.
2263 */
2264
2265 if (printlevel >= 2)
2266 printf("Total Iterations: %d \n",iter);
2267 /*
2268 * Now, go ahead and return.
2269 */
2270
2271 return(retcode);
2272
2273 }
2274
2275
2276