1 /*
2 * Copyright (c) 1994 by Lieven Vandenberghe and Stephen Boyd.
3 * Permission to use, copy, modify, and distribute this software for
4 * any purpose without fee is hereby granted, provided that this entire
5 * notice is included in all copies of any software which is or includes
6 * a copy or modification of this software and in all copies of the
7 * supporting documentation for such software.
8 * This software is being provided "as is", without any express or
9 * implied warranty. In particular, the authors do not make any
10 * representation or warranty of any kind concerning the merchantability
11 * of this software or its fitness for any particular purpose.
12 */
13 /*--------------------------------------------------------------------------*/
14 #include <stdio.h>
15 #include <math.h>
16 #include <string.h>
17 #include "spd.h"
18 #include "sciprint.h"
19 #include "localization.h"
20 #include "configvariable_interface.h"
21 #include "numericconstants_interface.h"
22 /*--------------------------------------------------------------------------*/
23 /* BLAS 1 */
24 extern double F2C(dnrm2)( );
25 extern double F2C(ddot)( );
26 extern void F2C(dcopy)( );
27 extern void F2C(daxpy)( );
28 extern void F2C(dscal)( );
29
30 /* BLAS 2 */
31 extern void F2C(dgemv)( );
32 extern void F2C(dspmv)( );
33
34 /* BLAS 3 */
35 extern void F2C(dgemm)( );
36
37 /* LAPACK */
38 extern void F2C(dgels)( );
39 extern void F2C(dspgst)( );
40 extern void F2C(dspev)( );
41 extern void F2C(dspgv)( );
42 extern void F2C(dtrcon)( );
43 /*--------------------------------------------------------------------------*/
44 /*
45 * if itype = 1, computes C = B*A*B', otherwise, computes C = B'*A*B
46 * A and B are nxn with A symmetric.
47 *
48 * Arguments:
49 * @param itype = 1: compute C = B*A*B'
50 * = any other integer: computes C = B'*A*B
51 * @param n dimension of A and B
52 * @param AP (input) double array of size n*(n+1)2;
53 * the lower triangle of A in packed storage
54 * @param B (input) double array of size n*n;
55 * @param CP (output) double array of size n*(n+1)/2;
56 * the lower triangle of C in packed storage
57 * @param temp: n-array, workspace
58 */
59 /*--------------------------------------------------------------------------*/
cngrncb(int itype,int n,double * AP,double * B,double * CP,double * temp)60 static void cngrncb(int itype, int n, double *AP, double *B, double *CP, double *temp)
61 {
62
63 int j, pos, lngth = n * (n + 1) / 2;
64 int int1 = 1;
65 double dbl0 = 0.0, dbl1 = 1.0;
66
67 /* C := 0 */
68 F2C(dscal)(&lngth, &dbl0, CP, &int1);
69
70 if (itype == 1)
71 {
72
73 for (j = 0, pos = 0; j < n; pos += n - j, j++)
74 {
75
76 /* temp = A*B(j,:)' */
77 F2C(dspmv)("L", &n, &dbl1, AP, B + j, &n, &dbl0, temp, &int1);
78
79 /* C(j:n,j) = B(j:n,:)*temp */
80 lngth = n - j;
81 F2C(dgemv)("N", &lngth, &n, &dbl1, B + j, &n, temp, &int1, &dbl0,
82 CP + pos, &int1);
83
84 }
85
86 }
87 else
88 {
89
90 for (j = 0, pos = 0; j < n; pos += n - j, j++)
91 {
92
93 /* temp = A*B(:,j) */
94 F2C(dspmv)("L", &n, &dbl1, AP, B + j * n, &int1, &dbl0, temp, &int1);
95
96 /* C(j:n,j) = B(:,j:n)'*temp */
97 lngth = n - j;
98 F2C(dgemv)("T", &n, &lngth, &dbl1, B + j * n, &n, temp, &int1, &dbl0,
99 CP + pos, &int1);
100
101 }
102 }
103
104 }
105 /*--------------------------------------------------------------------------*/
inprd(double * X,double * Z,int L,int * blck_szs)106 static double inprd(double *X, double *Z, int L, int *blck_szs)
107 /*
108 * Computes Tr X*Z
109 *
110 * Arguments:
111 * X,Z: block diagonal matrices with L blocks X^0, ..., X^{L-1},
112 * and Z^0, ..., Z^{L-1}. X^j and Z^j have size
113 * blck_szs[j] times blck_szs[j]. Every block is stored
114 * using packed storage of the lower triangle.
115 * L: number of blocks
116 * blck_szs: int vector of length L
117 * blck_szs[i], i=0,...,L-1 is the size of block i
118 *
119 */
120
121 {
122 double result;
123 int i, j, k, lngth, pos, sz, int1 = 1;
124
125 /* sz = length of Z and X */
126 for (i = 0, sz = 0; i < L; i++)
127 {
128 sz += (blck_szs[i] * (blck_szs[i] + 1)) / 2;
129 }
130
131 /* result = Tr X Z + contributions of diagonal elements */
132 result = 2.0 * F2C(ddot)(&sz, X, &int1, Z, &int1);
133
134 /* correct for diagonal elements
135 * loop over blocks, j=0,...,L-1 */
136 for (j = 0, pos = 0; j < L; j++)
137
138 /* loop over columns, k=0,...,blck_szs[j]-1
139 * pos is position of (k,k) element of block j
140 * lngth is length of column k */
141 for (k = 0, lngth = blck_szs[j]; k < blck_szs[j]; pos += lngth,
142 lngth -= 1, k++)
143
144 /* subtract Z^j_{kk}*X^j_{kk} from result */
145 {
146 result -= Z[pos] * X[pos];
147 }
148
149 return result;
150 }
151 /*--------------------------------------------------------------------------*/
C2F(spf)152 int C2F(spf)(
153 int *m, /* no of variables */
154 int *L, /* no of blocks in F */
155 double *F, /* F_i's in packed storage */
156 int *blck_szs, /* L-vector, dimensions of diagonal blocks */
157 double *c, /* m-vector */
158 double *x, /* m-vector */
159 double *Z, /* block diagonal matrix in packed storage */
160 double *ul, /* ul[0] = pr. obj, ul[1] = du. obj */
161 double *nu, /* >= 1.0 */
162 double *abstol, /* absolute accuracy */
163 double *reltol, /* relative accuracy */
164 double *tv, /* target value */
165 int *iters, /* on entry: the maximum number of iterations,
166 * on exit: the number of iterations taken */
167 double *work, /* work array */
168 int *lwork, /* size of work */
169 int *iwork, /* work array of m integers */
170 int *info /* status on termination */
171 )
172 {
173 return (sp(*m, *L, F, blck_szs, c, x, Z, ul, *nu, *abstol, *reltol, *tv, iters, work,
174 *lwork, iwork, info));
175 }
176 /*--------------------------------------------------------------------------*/
sp(int m,int L,double * F,int * blck_szs,double * c,double * x,double * Z,double * ul,double nu,double abstol,double reltol,double tv,int * iters,double * work,int lwork,int * iwork,int * info)177 int sp(
178 int m, /* no of variables */
179 int L, /* no of blocks in F */
180 double *F, /* F_i's in packed storage */
181 int *blck_szs, /* L-vector, dimensions of diagonal blocks */
182 double *c, /* m-vector */
183 double *x, /* m-vector */
184 double *Z, /* block diagonal matrix in packed storage */
185 double *ul, /* ul[0] = pr. obj, ul[1] = du. obj */
186 double nu, /* >= 1.0 */
187 double abstol, /* absolute accuracy */
188 double reltol, /* relative accuracy */
189 double tv, /* target value */
190 int *iters, /* on entry: the maximum number of iterations,
191 * on exit: the number of iterations taken */
192 double *work, /* work array */
193 int lwork, /* size of work */
194 int *iwork, /* work array of m integers */
195 int *info /* status on termination */
196 )
197 /*
198 * Solves semidefinite program
199 *
200 * minimize c'*x
201 * subject to F_0 + x_1*F_1 + ... + x_m*F_m >= 0
202 *
203 * and its dual
204 *
205 * maximize -Tr F_0*Z
206 * subject to Z >= 0
207 * Tr F_i*Z = c_i, i=1,...,m
208 *
209 *
210 * Convergence criteria:
211 * (1) maxiters is exceeded
212 * (2) duality gap is less than abstol
213 * (3) primal and dual objective are both positive and
214 * duality gap is less than reltol * dual objective
215 * or primal and dual objective are both negative and
216 * duality gap is less than reltol * minus the primal objective
217 * (4) reltol is negative and primal objective is less than tv
218 * (5) reltol is negative and dual objective is greater than tv
219 *
220 * Arguments:
221 * - m: number of variables x_i. m >= 1.
222 * - L: number of diagonal blocks in F_i. L >= 1.
223 * - F: the block diagonal matrices F_i, i=0,...,m.
224 * it is assumed that the matrices F_i are linearly
225 * independent.
226 * let F_i^j, i=0,..,m, j=0,...,L-1 denote the jth
227 * diagonal block of F_i,
228 * the array F contains F_0^0, ..., F_0^{L-1}, F_1^0, ...,
229 * F_1^{L-1}, ..., F_m^0, ..., F_m^{L-1}, in this order,
230 * using packed storage for the lower triangular part of
231 * F_i^j.
232 * - blck_szs: an int L-vector. blck_szs[j], j=0,....L-1 gives the
233 * size of block j, ie, F_i^j has size blck_szs[j]
234 * times blck_szs[j].
235 * - c: m-vector, primal objective.
236 * - x: m-vector. On entry, a strictly primal feasible point.
237 * On exit, the last iterate for x.
238 * - Z: block diagonal matrix with L blocks Z^0, ..., Z^{L-1}.
239 * Z^j has size blck_szs[j] times blck_szs[j].
240 * Every block is stored using packed storage of the lower
241 * triangular part.
242 * On entry, a strictly dual feasible point. On exit, the
243 * last dual iterate.
244 * - ul: two-vector. On exit, ul[0] is the primal objective value
245 * c'*x; ul[1] is the dual objective value -Tr F_0*Z.
246 * - nu: >= 1.0. Controls the rate of convergence.
247 * - abstol: absolute tolerance, >= MINABSTOL.
248 * - reltol: relative tolerance. Has a special meaning when negative.
249 * - tv: target value, only referenced if reltol < 0.
250 * - iters: on entry: maximum number of iterations >= 0,
251 * on exit: the number of iterations taken.
252 * - work: work array of size lwork.
253 * - lwork: size of work, must be at least:
254 * (m+2)*sz + up_sz + 2*n + ltemp, with
255 * ltemp = max( m+sz*nb, 3max_n + max_n*(max_n+1), 3*m )
256 * (sz: space needed to store one matrix F_i in packed
257 * storage, ie,
258 * sum_{j=0}^{L-1} blck_szs[j]*(blck_szs[j]+1)/2;
259 * up_sz: space needed to store one matrix F_i in
260 * unpacked storage, ie,
261 * sum_{j=0}^{L-1} blck_szs[j]*blck_szs[j];
262 * max_n: max block size;
263 * n: sum of the block sizes.
264 * nb >= 1, for best performance, nb should be at least
265 * equal to the optimal block size for dgels.
266 * - iwork: work array of m integers
267 * - info: returns 1 if maxiters exceeded, 2 if absolute accuracy
268 * is reached, 3 if relative accuracy is reached,
269 * 4 if target value is reached, 5 if target value is
270 * not achievable;
271 * negative values indicate errors: -i means argument i
272 * has an illegal value, -18 stands for all other errors.
273 *
274 *
275 * Returns 0 for normal exit, 1 if an error occurred.
276 *
277 */
278
279
280 {
281 int i, j, k, n, sz, up_sz, max_n, lngth, pos, pos2, pos3, pos4, ltemp,
282 maxiters, info2, minlwork;
283 double q, *rhs, *Fsc, *R, *X, rho, *dx, *sigx, *sigz, *dZ, *temp, scal,
284 scal2, XdZ, ZdX, alphax, alphaz, lambda_ls, gradx, hessx,
285 gradz, hessz, dalphax, dalphaz, gap, newgap = 0.0, newu = 0.0,
286 newl = 0.0, maxpossigx, minnegsigx, maxpossigz, minnegsigz, nrmc,
287 nrmx, nrmz, nrmmax, rcond;
288 int int2 = 2, int1 = 1;
289 double dbl1 = 1.0, dbl0 = 0.0, sqrt2 = sqrt(2.0);
290 double dbl_epsilon;
291
292 if (m < 1)
293 {
294 if ( getWarningMode() )
295 {
296 sciprint(_("m must be at least one.\n"));
297 }
298 *info = -1;
299 return 1;
300 }
301 if (L < 1)
302 {
303 if ( getWarningMode() )
304 {
305 sciprint(_("L must be at least one.\n"));
306 }
307 *info = -2;
308 return 1;
309 }
310 for (i = 0; i < L; i++) if (blck_szs[i] < 1)
311 {
312 if ( getWarningMode() )
313 {
314 sciprint(_("blck_szs[%d] must be at least one.\n"), i);
315 }
316 *info = -4;
317 return 1;
318 }
319 if (nu < 1.0)
320 {
321 if ( getWarningMode() )
322 {
323 sciprint(_("nu must be at least 1.0.\n"));
324 }
325 *info = -9;
326 return 1;
327 }
328
329
330 /*
331 * calculate dimensions:
332 * n: total size of semidefinite program
333 * sz: length of one block-diagonal matrix in packed storage
334 * up_sz: length of one block-diagonal matrix in unpacked storage
335 * max_n: size of biggest block
336 */
337
338 for (i = 0, n = 0, sz = 0, up_sz = 0, max_n = 0; i < L; i++)
339 {
340 n += blck_szs[i];
341 sz += blck_szs[i] * (blck_szs[i] + 1) / 2;
342 up_sz += blck_szs[i] * blck_szs[i];
343 max_n = Max(max_n, blck_szs[i]);
344 }
345 if (m > sz)
346 {
347 sciprint(_("Matrices Fi, i=1,...,m are linearly dependent.\n"));
348 *info = -3;
349 return 1;
350 }
351
352 q = (double)n + nu * sqrt((double)n);
353
354
355 /*
356 * check if Tr Fi*Z = c_i, i=1,...,m
357 */
358
359 nrmc = F2C(dnrm2)(&m, c, &int1);
360 for (i = 0; i < m; i++)
361 if (fabs(inprd(F + (i + 1)*sz, Z, L, blck_szs) - c[i]) > nrmc * TOLC)
362 {
363 if ( getWarningMode() )
364 {
365 sciprint(_("Z0 does not satisfy equality conditions for dual feasibility.\n"));
366 }
367
368 *info = -7;
369 return 1;
370 }
371
372
373 /*
374 * organize workspace
375 *
376 * work: (m+2)*sz + up_sz + 2*n + ltemp
377 * minimum ltemp: the maximum of
378 * m+sz*nb, 3*max_n + max_n*(max_n+1), and 3*m
379 * (nb is at least one)
380 *
381 * for dgels: m + sz*nb, nb at least 1
382 * for dspev("N"): 3*max_n + max_n*(max_n+1)
383 * for dspgv("N"): 3*max_n + max_n*(max_n+1)
384 * for dspgv("V"): 3*max_n + max_n*(max_n+1)/2
385 * for cngrncb: max_n
386 * for dtrcon: 3*m
387 *
388 * rhs (sz): work[0 ... sz-1]
389 * Fsc (m*sz): work[sz ... (m+1)*sz-1]
390 * R (up_sz): work[(m+1)*sz ... (m+1)*sz+up_sz-1]
391 * X (sz): work[(m+1)*sz+up_sz ... (m+2)*sz+up_sz-1]
392 * sigx (n): work[(m+2)*sz+up_sz ... (m+2)*sz+up_sz+n-1]
393 * sigz (n): work[(m+2)*sz+up_sz+n ... (m+2)*sz+up_sz+2*n-1]
394 * temp (remainder): work[(m+2)*sz+up_sz+2*n ... lwork-1]
395 */
396
397 /* check lwork */
398 minlwork = (m + 2) * sz + up_sz + 2 * n +
399 Max( Max( m + sz, 3 * max_n + max_n * (max_n + 1) ), 3 * m );
400 if (lwork < minlwork)
401 {
402 if ( getWarningMode() )
403 {
404 sciprint(_("Work space is too small. Need at least %d*sizeof(double).\n"), minlwork);
405 }
406 *info = -15;
407 return 1;
408 }
409
410 rhs = work; /* rhs for ls problem */
411 dx = work; /* solution of ls system; overlaps with rhs */
412 Fsc = rhs + sz; /* scaled matrices */
413 dZ = rhs + sz; /* overlaps with first column of Fsc */
414 R = Fsc + m * sz; /* eigenvectors of Z*F */
415 X = R + up_sz; /* F(x) */
416 sigx = X + sz; /* generalized eigenvalues of (dX,X) */
417 sigz = sigx + n; /* generalized eigenvalues of (dZ,Z) */
418 temp = sigz + n;
419 ltemp = lwork - (m + 2) * sz - up_sz - 2 * n;
420
421
422 maxiters = (*iters >= 0) ? *iters : MAXITERS;
423 for (*iters = 0; *iters <= maxiters; (*iters)++)
424 {
425
426
427 /* compute F(x) = F_0 + x_1*F_1 + ... + x_m*F_m, store in X */
428 F2C(dcopy)(&sz, F, &int1, X, &int1);
429 F2C(dgemv)("N", &sz, &m, &dbl1, F + sz, &sz, x, &int1, &dbl1, X, &int1);
430
431
432 /*
433 * compute generalized eigendecomp Z*F*x = lambda*x
434 * loop over blocks, i=0,...,L-1
435 * pos: position of (0,0) element of block i in packed storage
436 * pos2: position of (0,0) element of block i in unpacked
437 * storage
438 * pos3: position of first eigenvalue of block i in sigx
439 */
440
441 for (i = 0, pos = 0, pos2 = 0, pos3 = 0, gap = 0.0; i < L;
442 pos += blck_szs[i] * (blck_szs[i] + 1) / 2,
443 pos2 += blck_szs[i] * blck_szs[i],
444 pos3 += blck_szs[i], i++)
445 {
446
447 lngth = blck_szs[i] * (blck_szs[i] + 1) / 2;
448
449 /* copy block i of Z in temp (need max_n*(max_n+1)/2) */
450 F2C(dcopy)(&lngth, Z + pos, &int1, temp, &int1);
451
452 /* generalized eigenvalue decomposition Z*F*x = lambda*x
453 * - eigenvectors V are normalized s.t. V^T*F*V = I
454 * - store block i of V in R+pos2
455 * - store eigenvalues of block i in sigx+pos3
456 * - dspgv replaces X+pos by cholesky factor L of ith
457 * block of F (F = L*L^T)
458 * use temp+lngth as workspace (need at least 3*max_n) */
459 F2C(dspgv)(&int2, "V", "L", blck_szs + i, temp, X + pos, sigx + pos3,
460 R + pos2, blck_szs + i, temp + lngth, &info2);
461 if (info2)
462 {
463 if ( getWarningMode() )
464 {
465 sciprint(_("Error in dspgv, info = %d.\n"), info2);
466 }
467 if (*iters == 0 && info2 > blck_szs[i])
468 {
469 if ( getWarningMode() )
470 {
471 sciprint( _("x0 is not strictly primal feasible.\n"));
472 }
473 *info = -6;
474 }
475 else
476 {
477 *info = -18;
478 }
479 return 1;
480 }
481
482 /* - replace sigx+pos3 by lambda^(1/2)
483 * - normalize block i of V (stored in R+pos2) s.t.
484 * V^T*F*V = Lambda^(1/2) */
485 for (k = 0; k < blck_szs[i]; k++)
486 {
487 scal = sigx[pos3 + k];
488 if (scal < 0.0)
489 {
490 if (*iters == 0)
491 {
492 if ( getWarningMode() )
493 {
494 sciprint(_("Z0 is not positive definite.\n"));
495 }
496 *info = 7;
497 }
498 else
499 {
500 if ( getWarningMode() )
501 {
502 sciprint(_("F(x)*Z has a negative eigenvalue.\n"));
503 }
504 *info = -18;
505 }
506 return 1;
507 }
508 gap += scal; /* duality gap is sum of eigenvalues of ZF */
509 scal2 = sqrt(scal);
510 scal = sqrt(scal2);
511 sigx[pos3 + k] = scal2;
512 F2C(dscal)(blck_szs + i, &scal, R + pos2 + k * blck_szs[i], &int1);
513 }
514
515 }
516
517
518 /*
519 * check convergence
520 */
521
522 ul[1] = -inprd(F, Z, L, blck_szs); /* -Tr F_0 Z */
523 ul[0] = F2C(ddot)(&m, c, &int1, x, &int1); /* c^T x */
524 if (*iters == 0)
525 {
526 if ( getWarningMode() )
527 {
528 sciprint(_("\n primal obj. dual obj. dual. gap \n"));
529 }
530
531 }
532 if ( getWarningMode() )
533 {
534 sciprint("% 13.2e % 12.2e %10.2e\n", ul[0], ul[1], gap);
535 }
536 if (gap <= Max(abstol, MINABSTOL))
537 {
538 *info = 2;
539 }
540 else if ( (ul[1] > 0.0 && gap <= reltol * ul[1]) ||
541 (ul[0] < 0.0 && gap <= reltol * (-ul[0])) )
542 {
543 *info = 3;
544 }
545 else if ( reltol < 0.0 && ul[0] <= tv )
546 {
547 *info = 4;
548 }
549 else if ( reltol < 0.0 && ul[1] >= tv )
550 {
551 *info = 5;
552 }
553 else if ( *iters == maxiters )
554 {
555 *info = 1;
556 }
557 else
558 {
559 *info = 0;
560 }
561 if (*info)
562 {
563 return 0;
564 }
565
566
567
568 /*
569 * compute scaled matrices F
570 */
571
572 for (j = 0, pos = 0; j < m; j++) for (i = 0, pos2 = 0; i < L;
573 pos += blck_szs[i] * (blck_szs[i] + 1) / 2,
574 pos2 += blck_szs[i] * blck_szs[i], i++)
575 {
576
577 /* compute R' * Fj(i) * R, store in Fsc+pos */
578 cngrncb(2, blck_szs[i], F + sz + pos, R + pos2, Fsc + pos, temp);
579
580 /* correct diagonal elements */
581 for (k = 0, pos4 = pos; k < blck_szs[i]; pos4 += blck_szs[i] - k, k++)
582 {
583 Fsc[pos4] /= sqrt2;
584 }
585
586 }
587
588
589 /*
590 * form rhs = Lambda^(-1/2) - (q/gap) * Lambda^(1/2)
591 */
592
593 F2C(dscal)(&sz, &dbl0, rhs, &int1); /* rhs := 0 */
594 rho = -q / gap;
595 for (i = 0, pos = 0, pos3 = 0; i < L;
596 pos += blck_szs[i] * (blck_szs[i] + 1) / 2,
597 pos3 += blck_szs[i], i++)
598 for (k = 0, pos4 = pos; k < blck_szs[i]; pos4 += blck_szs[i] - k, k++)
599 {
600 scal = sigx[pos3 + k];
601 rhs[pos4] = (1.0 / scal + rho * scal) / sqrt2;
602 }
603
604
605 /*
606 * solve least-squares problem; need workspace of size m + nb*sz
607 * - rhs is overwritten by dx
608 * - in first iteration, estimate condition number of Fsc
609 */
610
611 F2C(dgels)("N", &sz, &m, &int1, Fsc, &sz, rhs, &sz, temp, <emp,
612 &info2);
613 if (info2)
614 {
615 if ( getWarningMode() )
616 {
617 sciprint(_("Error in dgels, info = %d.\n"), info2);
618 }
619 *info = -18;
620 return 1;
621 }
622
623 if (*iters == 0)
624 {
625
626 /* estimate the condition number in 1-norm of the R-factor of
627 * the qr-decomposition of Fsc (is stored in Fsc)
628 * need work space of size 3*m */
629 F2C(dtrcon)("1", "U", "N", &m, Fsc, &sz, &rcond, temp, iwork,
630 &info2);
631 if (info2 < 0)
632 {
633 if ( getWarningMode() )
634 {
635 sciprint(_("Error in dtrcon, info = %d.\n"), info2);
636 }
637 *info = -18;
638 return 1;
639 }
640 if (rcond < MINRCOND)
641 {
642 if ( getWarningMode() )
643 {
644 sciprint(_("The matrices F_i, i=1,...,m are linearly dependent (or the initial points are very badly conditioned).\n"));
645 }
646 *info = -3;
647 return 1;
648 }
649
650 }
651
652
653
654 /*
655 * - compute dZ =
656 * R*((q/gap)*Lambda^(1/2) - Lambda^(-1/2) + R^T*dF*R )*R^T
657 * - compute generalized eigenvalues of (dF, F), store in sigx
658 * - compute generalized eigenvalues of (dZ, Z), store in sigz
659 *
660 * loop over blocks i=0,...,L-1
661 * pos: position of (0,0) element of block i in packed storage
662 * pos2: position of (0,0) element of block i in unpacked storage
663 * pos3: position of first eigenvalue of in sigx and sigz
664 */
665
666 for (i = 0, pos = 0, pos2 = 0, pos3 = 0; i < L;
667 pos += blck_szs[i] * (blck_szs[i] + 1) / 2,
668 pos2 += blck_szs[i] * blck_szs[i],
669 pos3 += blck_szs[i], i++)
670 {
671
672 lngth = blck_szs[i] * (blck_szs[i] + 1) / 2;
673
674 /* compute ith block of dF = \sum \delta x_i F_i,
675 * store in temp */
676 F2C(dgemv)("N", &lngth, &m, &dbl1, F + sz + pos, &sz, dx, &int1,
677 &dbl0, temp, &int1);
678
679 /* scale dF as R'*dF*R, store in temp + lngth */
680 cngrncb(2, blck_szs[i], temp, R + pos2, temp + lngth, temp + 2 * lngth);
681
682 /* add (q/gap)*Lambda^(1/2) - Lambda^(-1/2) */
683 for (k = 0, pos4 = lngth; k < blck_szs[i]; pos4 += blck_szs[i] - k, k++)
684 {
685 temp[pos4] -= rho * sigx[pos3 + k] + 1.0 / sigx[pos3 + k];
686 }
687
688 /* replace dF in temp by L^{-1}*dF*L^{-T},
689 * (L: cholesky factor of F, stored in X)
690 * and compute eigenvalues of L^{-1}*dF*L^{-T} */
691 F2C(dspgst)(&int1, "L", blck_szs + i, temp, X + pos, &info2);
692 if (info2)
693 {
694 if ( getWarningMode() )
695 {
696 sciprint(_("Error in dspst, info = %d.\n"), info2);
697 }
698
699 *info = -18;
700 return 1;
701 }
702 /* temp has to be of size max_n*(max_n+1)+3*max_n */
703 F2C(dspev)("N", "L", blck_szs + i, temp, sigx + pos3, NULL, &int1,
704 temp + 2 * lngth, &info2);
705 if (info2)
706 {
707 if ( getWarningMode() )
708 {
709 sciprint(_("Error in dspev, info = %d.\n"), info2);
710 }
711 *info = -18;
712 return 1;
713 }
714
715 /* dZ := R*((q/gap)*Lambda^(1/2) - Lambda^(-1/2) + R'*dF*R)*R' */
716 cngrncb(1, blck_szs[i], temp + lngth, R + pos2, dZ + pos,
717 temp + 2 * lngth);
718
719 /* copy ith block of dZ to temp */
720 F2C(dcopy)(&lngth, dZ + pos, &int1, temp, &int1);
721
722 /* copy ith block of Z to temp + lngth */
723 F2C(dcopy)(&lngth, Z + pos, &int1, temp + lngth, &int1);
724
725 /* sigz: generalized eigenvalues of (dZ,Z)
726 * required size of temp: 3*max_n + max_n*(max_n+1) */
727 F2C(dspgv)(&int1, "N", "L", blck_szs + i, temp, temp + lngth, sigz + pos3,
728 NULL, &int1, temp + 2 * lngth, &info2);
729 if (info2)
730 {
731 if ( getWarningMode() )
732 {
733 sciprint(_("Error in dspgv, info = %d.\n"), info2);
734 }
735 *info = -18;
736 return 1;
737 }
738
739 }
740
741
742 /*
743 * compute feasible rectangle for plane search
744 */
745
746 maxpossigx = 0.0;
747 minnegsigx = 0.0;
748 maxpossigz = 0.0;
749 minnegsigz = 0.0;
750 for (i = 0; i < n; i++)
751 {
752 if ( sigx[i] > maxpossigx )
753 {
754 maxpossigx = sigx[i]; /* max pos eigenvalue in sigx */
755 }
756 else if ( sigx[i] < minnegsigx )
757 {
758 minnegsigx = sigx[i]; /* min neg eigenvalue in sigx */
759 }
760 if ( sigz[i] > maxpossigz )
761 {
762 maxpossigz = sigz[i]; /* max pos eigenvalue in sigz */
763 }
764 else if ( sigz[i] < minnegsigz )
765 {
766 minnegsigz = sigz[i]; /* min neg eigenvalue in sigz */
767 }
768 }
769 nrmx = F2C(dnrm2)(&n, sigx, &int1); /* norm of scaled dx */
770 nrmz = F2C(dnrm2)(&n, sigz, &int1); /* norm of scaled dZ */
771 nrmmax = Max( nrmx, nrmz);
772
773 XdZ = inprd(F, dZ, L, blck_szs); /* Tr F0*dZ */
774 ZdX = F2C(ddot)(&m, c, &int1, dx, &int1); /* c^T*dx */
775
776
777 /*
778 * check corners of feasible rectangle
779 */
780
781 dbl_epsilon = nc_eps();
782 if (nrmx > SIGTOL * nrmmax)
783 if (ZdX < 0.0)
784 {
785 alphax = (minnegsigx < -dbl_epsilon) ? -1.0 / minnegsigx : 0.0;
786 }
787 else
788 {
789 alphax = (maxpossigx > dbl_epsilon) ? -1.0 / maxpossigx : 0.0;
790 }
791 else
792 {
793 alphax = 0.0;
794 }
795
796 if (nrmz > SIGTOL * nrmmax)
797 if (XdZ < 0.0)
798 {
799 alphaz = (minnegsigz < -dbl_epsilon) ? -1.0 / minnegsigz : 0.0;
800 }
801 else
802 {
803 alphaz = (maxpossigz > dbl_epsilon) ? -1.0 / maxpossigz : 0.0;
804 }
805 else
806 {
807 alphaz = 0.0;
808 }
809
810 newgap = gap + alphax * ZdX + alphaz * XdZ;
811 newu = ul[0] + alphax * ZdX;
812 newl = ul[1] - alphaz * XdZ;
813
814 if (newgap <= Max(abstol, MINABSTOL))
815 {
816 *info = 2;
817 }
818 else if ( (newl > 0.0 && newgap <= reltol * newl) ||
819 (newu < 0.0 && newgap <= -reltol * newu) )
820 {
821 *info = 3;
822 }
823 else if ( reltol < 0.0 && newu <= tv )
824 {
825 *info = 4;
826 }
827 else if ( reltol < 0.0 && newl >= tv )
828 {
829 *info = 5;
830 }
831 else if ( *iters == maxiters )
832 {
833 *info = 1;
834 }
835 else
836 {
837 *info = 0;
838 }
839
840 if (*info)
841 {
842 F2C(daxpy)(&m, &alphax, dx, &int1, x, &int1);
843 F2C(daxpy)(&sz, &alphaz, dZ, &int1, Z, &int1);
844 gap = newgap;
845 ul[0] = newu;
846 ul[1] = newl;
847 if ( getWarningMode() )
848 {
849 sciprint("% 13.2e % 12.2e %10.2e\n", ul[0], ul[1], gap);
850 }
851 (*iters)++;
852 return 0;
853 }
854
855
856 /*
857 * plane search
858 * minimize phi(alphax,alphaz) =
859 * q*log(dual_gap + alphax*c^T*dx + alphaz* Tr F_0 dZ)
860 * - sum log (1+alphax*sigx_i) - sum log (1+alphaz*sigz)
861 */
862
863 alphax = 0.0;
864 alphaz = 0.0;
865 lambda_ls = 1.0;
866
867 if (nrmx > SIGTOL * nrmmax)
868 if (nrmz > SIGTOL * nrmmax) /* compute primal and dual steps */
869 while ( lambda_ls > 1e-4 )
870 {
871
872 /* compute 1st and 2nd derivatives of phi */
873 rho = q / (gap + alphax * ZdX + alphaz * XdZ);
874 gradx = rho * ZdX;
875 hessx = 0.0;
876 gradz = rho * XdZ;
877 hessz = 0.0;
878 for (i = 0; i < n; i++)
879 {
880 gradx -= sigx[i] / (1.0 + alphax * sigx[i]);
881 hessx += SQR( sigx[i] / (1.0 + alphax * sigx[i]) );
882 gradz -= sigz[i] / (1.0 + alphaz * sigz[i]);
883 hessz += SQR( sigz[i] / (1.0 + alphaz * sigz[i]) );
884 }
885
886 /* newton step */
887 dalphax = -gradx / hessx;
888 dalphaz = -gradz / hessz;
889 lambda_ls = sqrt( SQR(gradx) / hessx + SQR(gradz) / hessz );
890 alphax += (lambda_ls > 0.25) ?
891 dalphax / (1.0 + lambda_ls) : dalphax;
892 alphaz += (lambda_ls > 0.25) ?
893 dalphaz / (1.0 + lambda_ls) : dalphaz;
894
895 }
896
897 else while ( lambda_ls > 1e-4 ) /* primal step only */
898 {
899
900 /* compute 1st and 2nd derivatives of phi */
901 rho = q / (gap + alphax * ZdX);
902 gradx = rho * ZdX;
903 hessx = 0.0;
904 for (i = 0; i < n; i++)
905 {
906 gradx -= sigx[i] / (1.0 + alphax * sigx[i]);
907 hessx += SQR( sigx[i] / (1.0 + alphax * sigx[i]) );
908 }
909
910 /* newton step */
911 dalphax = -gradx / hessx;
912 lambda_ls = fabs(gradx) / sqrt(hessx);
913 alphax += (lambda_ls > 0.25) ?
914 dalphax / (1.0 + lambda_ls) : dalphax;
915
916 }
917
918 else if (nrmz > SIGTOL * nrmmax) /* dual step only */
919 while ( lambda_ls > 1e-4 )
920 {
921
922 /* compute 1st and 2nd derivatives of phi */
923 rho = q / (gap + alphaz * XdZ);
924 gradz = rho * XdZ;
925 hessz = 0.0;
926 for (i = 0; i < n; i++)
927 {
928 gradz -= sigz[i] / (1.0 + alphaz * sigz[i]);
929 hessz += SQR( sigz[i] / (1.0 + alphaz * sigz[i]) );
930 }
931
932 /* newton step */
933 dalphaz = -gradz / hessz;
934 lambda_ls = fabs(gradz) / sqrt(hessz);
935 alphaz += (lambda_ls > 0.25) ?
936 dalphaz / (1.0 + lambda_ls) : dalphaz;
937 }
938
939
940
941 /* update x and Z */
942 F2C(daxpy)(&m, &alphax, dx, &int1, x, &int1);
943 F2C(daxpy)(&sz, &alphaz, dZ, &int1, Z, &int1);
944
945 }
946
947 return -1; /* should never happen */
948 }
949
950