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.
19
void cngrncb(itype,n,AP,B,CP,temp)
28
* if itype = 1, computes C = B*A*B', otherwise, computes C = B'*A*B
29
* A and B are nxn with A symmetric.
32
* - itype = 1: compute C = B*A*B'
33
* = any other integer: computes C = B'*A*B
34
* - n dimension of A and B
35
* - AP (input) double array of size n*(n+1)2;
36
* the lower triangle of A in packed storage
37
* - B (input) double array of size n*n;
38
* - CP (output) double array of size n*(n+1)/2;
39
* the lower triangle of C in packed storage
40
* - temp: n-array, workspace
45
int j, pos, lngth = n*(n+1)/2;
47
double dbl0=0.0, dbl1=1.0;
50
F2C(dscal)(&lngth, &dbl0, CP, &int1);
54
for (j=0, pos=0; j<n; pos+=n-j, j++){
56
/* temp = A*B(j,:)' */
57
F2C(dspmv)("L", &n, &dbl1, AP, B+j, &n, &dbl0, temp, &int1);
59
/* C(j:n,j) = B(j:n,:)*temp */
61
F2C(dgemv)("N", &lngth, &n, &dbl1, B+j, &n, temp, &int1, &dbl0,
68
for (j=0, pos=0; j<n; pos+=n-j, j++){
71
F2C(dspmv)("L", &n, &dbl1, AP, B+j*n, &int1, &dbl0, temp, &int1);
73
/* C(j:n,j) = B(:,j:n)'*temp */
75
F2C(dgemv)("T", &n, &lngth, &dbl1, B+j*n, &n, temp, &int1, &dbl0,
84
double inprd(X,Z,L,blck_szs)
94
* X,Z: block diagonal matrices with L blocks X^0, ..., X^{L-1},
95
* and Z^0, ..., Z^{L-1}. X^j and Z^j have size
96
* blck_szs[j] times blck_szs[j]. Every block is stored
97
* using packed storage of the lower triangle.
99
* blck_szs: integer vector of length L
100
* blck_szs[i], i=0,...,L-1 is the size of block i
106
int i, j, k, lngth, pos, sz, int1=1;
108
/* sz = length of Z and X */
109
for (i=0, sz=0; i<L; i++) sz += (blck_szs[i]*(blck_szs[i]+1))/2;
111
/* result = Tr X Z + contributions of diagonal elements */
112
result = 2.0*F2C(ddot)(&sz, X, &int1, Z, &int1);
114
/* correct for diagonal elements
115
* loop over blocks, j=0,...,L-1 */
116
for (j=0, pos=0; j<L; j++)
118
/* loop over columns, k=0,...,blck_szs[j]-1
119
* pos is position of (k,k) element of block j
120
* lngth is length of column k */
121
for (k=0, lngth=blck_szs[j]; k<blck_szs[j]; pos+=lngth,
124
/* subtract Z^j_{kk}*X^j_{kk} from result */
125
result -= Z[pos]*X[pos];
130
int C2F(spf)(m,L,F,blck_szs,c,x,Z,ul,nu,abstol,reltol,tv,iters,
131
work,lwork,iwork,info)
132
int *m; /* no of variables */
133
int *L; /* no of blocks in F */
134
double *F; /* F_i's in packed storage */
135
int *blck_szs; /* L-vector, dimensions of diagonal blocks */
136
double *c; /* m-vector */
137
double *x; /* m-vector */
138
double *Z; /* block diagonal matrix in packed storage */
139
double *ul; /* ul[0] = pr. obj, ul[1] = du. obj */
140
double *nu; /* >= 1.0 */
141
double *abstol; /* absolute accuracy */
142
double *reltol; /* relative accuracy */
143
double *tv; /* target value */
144
int *iters; /* on entry: the maximum number of iterations,
145
* on exit: the number of iterations taken */
146
double *work; /* work array */
147
int *lwork; /* size of work */
148
int *iwork; /* work array of m integers */
149
int *info; /* status on termination */
151
return(sp(*m,*L,F,blck_szs,c,x,Z,ul,*nu,*abstol,*reltol,*tv,iters, work,
155
int sp(m,L,F,blck_szs,c,x,Z,ul,nu,abstol,reltol,tv,iters,work,
158
int m; /* no of variables */
159
int L; /* no of blocks in F */
160
double *F; /* F_i's in packed storage */
161
int *blck_szs; /* L-vector, dimensions of diagonal blocks */
162
double *c; /* m-vector */
163
double *x; /* m-vector */
164
double *Z; /* block diagonal matrix in packed storage */
165
double *ul; /* ul[0] = pr. obj, ul[1] = du. obj */
166
double nu; /* >= 1.0 */
167
double abstol; /* absolute accuracy */
168
double reltol; /* relative accuracy */
169
double tv; /* target value */
170
int *iters; /* on entry: the maximum number of iterations,
171
* on exit: the number of iterations taken */
172
double *work; /* work array */
173
int lwork; /* size of work */
174
int *iwork; /* work array of m integers */
175
int *info; /* status on termination */
178
* Solves semidefinite program
181
* subject to F_0 + x_1*F_1 + ... + x_m*F_m >= 0
187
* Tr F_i*Z = c_i, i=1,...,m
190
* Convergence criteria:
191
* (1) maxiters is exceeded
192
* (2) duality gap is less than abstol
193
* (3) primal and dual objective are both positive and
194
* duality gap is less than reltol * dual objective
195
* or primal and dual objective are both negative and
196
* duality gap is less than reltol * minus the primal objective
197
* (4) reltol is negative and primal objective is less than tv
198
* (5) reltol is negative and dual objective is greater than tv
201
* - m: number of variables x_i. m >= 1.
202
* - L: number of diagonal blocks in F_i. L >= 1.
203
* - F: the block diagonal matrices F_i, i=0,...,m.
204
* it is assumed that the matrices F_i are linearly
206
* let F_i^j, i=0,..,m, j=0,...,L-1 denote the jth
207
* diagonal block of F_i,
208
* the array F contains F_0^0, ..., F_0^{L-1}, F_1^0, ...,
209
* F_1^{L-1}, ..., F_m^0, ..., F_m^{L-1}, in this order,
210
* using packed storage for the lower triangular part of
212
* - blck_szs: an integer L-vector. blck_szs[j], j=0,....L-1 gives the
213
* size of block j, ie, F_i^j has size blck_szs[j]
215
* - c: m-vector, primal objective.
216
* - x: m-vector. On entry, a strictly primal feasible point.
217
* On exit, the last iterate for x.
218
* - Z: block diagonal matrix with L blocks Z^0, ..., Z^{L-1}.
219
* Z^j has size blck_szs[j] times blck_szs[j].
220
* Every block is stored using packed storage of the lower
222
* On entry, a strictly dual feasible point. On exit, the
224
* - ul: two-vector. On exit, ul[0] is the primal objective value
225
* c'*x; ul[1] is the dual objective value -Tr F_0*Z.
226
* - nu: >= 1.0. Controls the rate of convergence.
227
* - abstol: absolute tolerance, >= MINABSTOL.
228
* - reltol: relative tolerance. Has a special meaning when negative.
229
* - tv: target value, only referenced if reltol < 0.
230
* - iters: on entry: maximum number of iterations >= 0,
231
* on exit: the number of iterations taken.
232
* - work: work array of size lwork.
233
* - lwork: size of work, must be at least:
234
* (m+2)*sz + up_sz + 2*n + ltemp, with
235
* ltemp = max( m+sz*nb, 3max_n + max_n*(max_n+1), 3*m )
236
* (sz: space needed to store one matrix F_i in packed
238
* sum_{j=0}^{L-1} blck_szs[j]*(blck_szs[j]+1)/2;
239
* up_sz: space needed to store one matrix F_i in
240
* unpacked storage, ie,
241
* sum_{j=0}^{L-1} blck_szs[j]*blck_szs[j];
242
* max_n: max block size;
243
* n: sum of the block sizes.
244
* nb >= 1, for best performance, nb should be at least
245
* equal to the optimal block size for dgels.
246
* - iwork: work array of m integers
247
* - info: returns 1 if maxiters exceeded, 2 if absolute accuracy
248
* is reached, 3 if relative accuracy is reached,
249
* 4 if target value is reached, 5 if target value is
251
* negative values indicate errors: -i means argument i
252
* has an illegal value, -18 stands for all other errors.
255
* Returns 0 for normal exit, 1 if an error occurred.
261
int i, j, k, n, sz, up_sz, max_n, lngth, pos, pos2, pos3, pos4, ltemp,
262
maxiters, info2, minlwork;
263
double q, *rhs, *Fsc, *R, *X, rho, *dx, *sigx, *sigz, *dZ, *temp, scal,
264
scal2, XdZ, ZdX, alphax, alphaz, lambda_ls, gradx, hessx,
265
gradz, hessz, dalphax, dalphaz, gap, newgap=0.0, newu=0.0,
266
newl=0.0, maxpossigx, minnegsigx, maxpossigz, minnegsigz, nrmc,
267
nrmx, nrmz, nrmmax, rcond;
269
double dbl1=1.0, dbl0=0.0, sqrt2=sqrt(2.0);
274
sprintf(str, "m must be at least one. \n");
280
sprintf(str, "L must be at least one. \n");
285
for (i=0; i<L; i++) if (blck_szs[i] < 1){
286
sprintf(str, "blck_szs[%d] must be at least one.\n", i);
292
sprintf(str, "nu must be at least 1.0. \n");
300
* calculate dimensions:
301
* n: total size of semidefinite program
302
* sz: length of one block-diagonal matrix in packed storage
303
* up_sz: length of one block-diagonal matrix in unpacked storage
304
* max_n: size of biggest block
307
for (i=0, n=0, sz=0, up_sz=0, max_n=0; i<L; i++){
309
sz += blck_szs[i]*(blck_szs[i]+1)/2;
310
up_sz += blck_szs[i]*blck_szs[i];
311
max_n = MAX(max_n, blck_szs[i]);
314
sprintf(str, "The matrices Fi, i=1,...,m are linearly dependent.\n");
316
*info = -3; return 1;
319
q = (double)n + nu*sqrt((double)n);
323
* check if Tr Fi*Z = c_i, i=1,...,m
326
nrmc = F2C(dnrm2)(&m, c, &int1);
328
if (fabs(inprd(F+(i+1)*sz, Z, L, blck_szs) - c[i]) > nrmc*TOLC){
329
sprintf(str, "Z0 does not satisfy equality conditions\
330
for dual feasibility.\n");
340
* work: (m+2)*sz + up_sz + 2*n + ltemp
341
* minimum ltemp: the maximum of
342
* m+sz*nb, 3*max_n + max_n*(max_n+1), and 3*m
343
* (nb is at least one)
345
* for dgels: m + sz*nb, nb at least 1
346
* for dspev("N"): 3*max_n + max_n*(max_n+1)
347
* for dspgv("N"): 3*max_n + max_n*(max_n+1)
348
* for dspgv("V"): 3*max_n + max_n*(max_n+1)/2
352
* rhs (sz): work[0 ... sz-1]
353
* Fsc (m*sz): work[sz ... (m+1)*sz-1]
354
* R (up_sz): work[(m+1)*sz ... (m+1)*sz+up_sz-1]
355
* X (sz): work[(m+1)*sz+up_sz ... (m+2)*sz+up_sz-1]
356
* sigx (n): work[(m+2)*sz+up_sz ... (m+2)*sz+up_sz+n-1]
357
* sigz (n): work[(m+2)*sz+up_sz+n ... (m+2)*sz+up_sz+2*n-1]
358
* temp (remainder): work[(m+2)*sz+up_sz+2*n ... lwork-1]
362
minlwork = (m+2)*sz + up_sz + 2*n +
363
MAX( MAX( m+sz, 3*max_n + max_n*(max_n+1) ), 3*m );
364
if (lwork < minlwork){
365
sprintf(str, "Work space is too small. Need at least\
366
%d*sizeof(double).\n", minlwork);
372
rhs = work; /* rhs for ls problem */
373
dx = work; /* solution of ls system; overlaps with rhs */
374
Fsc = rhs + sz; /* scaled matrices */
375
dZ = rhs + sz; /* overlaps with first column of Fsc */
376
R = Fsc + m*sz; /* eigenvectors of Z*F */
377
X = R + up_sz; /* F(x) */
378
sigx = X + sz; /* generalized eigenvalues of (dX,X) */
379
sigz = sigx + n; /* generalized eigenvalues of (dZ,Z) */
381
ltemp = lwork - (m+2)*sz - up_sz - 2*n;
384
maxiters = (*iters >= 0) ? *iters : MAXITERS;
385
for (*iters=0; *iters <= maxiters; (*iters)++){
388
/* compute F(x) = F_0 + x_1*F_1 + ... + x_m*F_m, store in X */
389
F2C(dcopy)(&sz, F, &int1, X, &int1);
390
F2C(dgemv)("N", &sz, &m, &dbl1, F+sz, &sz, x, &int1, &dbl1, X, &int1);
394
* compute generalized eigendecomp Z*F*x = lambda*x
395
* loop over blocks, i=0,...,L-1
396
* pos: position of (0,0) element of block i in packed storage
397
* pos2: position of (0,0) element of block i in unpacked
399
* pos3: position of first eigenvalue of block i in sigx
402
for (i=0, pos=0, pos2=0, pos3=0, gap=0.0; i<L;
403
pos += blck_szs[i]*(blck_szs[i]+1)/2,
404
pos2 += blck_szs[i]*blck_szs[i],
405
pos3 += blck_szs[i], i++){
407
lngth = blck_szs[i]*(blck_szs[i]+1)/2;
409
/* copy block i of Z in temp (need max_n*(max_n+1)/2) */
410
F2C(dcopy)(&lngth, Z+pos, &int1, temp, &int1);
412
/* generalized eigenvalue decomposition Z*F*x = lambda*x
413
* - eigenvectors V are normalized s.t. V^T*F*V = I
414
* - store block i of V in R+pos2
415
* - store eigenvalues of block i in sigx+pos3
416
* - dspgv replaces X+pos by cholesky factor L of ith
417
* block of F (F = L*L^T)
418
* use temp+lngth as workspace (need at least 3*max_n) */
419
F2C(dspgv)(&int2, "V", "L", blck_szs+i, temp, X+pos, sigx+pos3,
420
R+pos2, blck_szs+i, temp+lngth, &info2);
422
sprintf(str,"Error in dspgv, info = %d.\n", info2);
424
if (*iters == 0 && info2 > blck_szs[i]){
425
sprintf(str, "x0 is not strictly primal feasible.\n");
432
/* - replace sigx+pos3 by lambda^(1/2)
433
* - normalize block i of V (stored in R+pos2) s.t.
434
* V^T*F*V = Lambda^(1/2) */
435
for (k=0; k<blck_szs[i]; k++){
439
sprintf(str, "Z0 is not positive definite.\n");
443
sprintf(str, "F(x)*Z has a negative eigenvalue.\n");
449
gap += scal; /* duality gap is sum of eigenvalues of ZF */
452
sigx[pos3+k] = scal2;
453
F2C(dscal)(blck_szs+i, &scal, R+pos2+k*blck_szs[i], &int1);
463
ul[1] = -inprd(F,Z,L,blck_szs); /* -Tr F_0 Z */
464
ul[0] = F2C(ddot)(&m, c, &int1, x, &int1); /* c^T x */
466
sprintf(str,"\n primal obj. dual obj. dual. gap \n");
469
sprintf(str,"% 13.2e % 12.2e %10.2e \n", ul[0], ul[1], gap);
471
if (gap <= MAX(abstol, MINABSTOL)) *info = 2;
472
else if ( (ul[1] > 0.0 && gap <= reltol*ul[1]) ||
473
(ul[0] < 0.0 && gap <= reltol*(-ul[0])) ) *info = 3;
474
else if ( reltol < 0.0 && ul[0] <= tv ) *info = 4;
475
else if ( reltol < 0.0 && ul[1] >= tv ) *info = 5;
476
else if ( *iters == maxiters ) *info = 1;
483
* compute scaled matrices F
486
for (j=0, pos=0; j<m; j++) for (i=0, pos2=0; i<L;
487
pos += blck_szs[i]*(blck_szs[i]+1)/2,
488
pos2 += blck_szs[i]*blck_szs[i], i++) {
490
/* compute R' * Fj(i) * R, store in Fsc+pos */
491
cngrncb(2, blck_szs[i], F+sz+pos, R+pos2, Fsc+pos, temp);
493
/* correct diagonal elements */
494
for (k=0, pos4=pos; k<blck_szs[i]; pos4 += blck_szs[i]-k, k++)
501
* form rhs = Lambda^(-1/2) - (q/gap) * Lambda^(1/2)
504
F2C(dscal)(&sz, &dbl0, rhs, &int1); /* rhs := 0 */
506
for (i=0, pos=0, pos3=0; i<L;
507
pos += blck_szs[i]*(blck_szs[i]+1)/2,
508
pos3 += blck_szs[i], i++)
509
for (k=0, pos4=pos; k<blck_szs[i]; pos4+=blck_szs[i]-k, k++){
511
rhs[pos4] = (1.0/scal + rho*scal)/sqrt2;
516
* solve least-squares problem; need workspace of size m + nb*sz
517
* - rhs is overwritten by dx
518
* - in first iteration, estimate condition number of Fsc
521
F2C(dgels)("N", &sz, &m, &int1, Fsc, &sz, rhs, &sz, temp, <emp,
524
sprintf(str,"Error in dgels, info = %d.\n", info2);
526
*info = -18; return 1;
531
/* estimate the condition number in 1-norm of the R-factor of
532
* the qr-decomposition of Fsc (is stored in Fsc)
533
* need work space of size 3*m */
534
F2C(dtrcon)("1", "U", "N", &m, Fsc, &sz, &rcond, temp, iwork,
537
sprintf(str,"Error in dtrcon, info = %d.\n", info2);
539
*info = -18; return 1;
541
if (rcond < MINRCOND) {
542
sprintf(str,"The matrices F_i, i=1,...,m are linearly\
543
dependent (or the initial points are very badly conditioned).\n");
545
*info = -3; return 1;
554
* R*((q/gap)*Lambda^(1/2) - Lambda^(-1/2) + R^T*dF*R )*R^T
555
* - compute generalized eigenvalues of (dF, F), store in sigx
556
* - compute generalized eigenvalues of (dZ, Z), store in sigz
558
* loop over blocks i=0,...,L-1
559
* pos: position of (0,0) element of block i in packed storage
560
* pos2: position of (0,0) element of block i in unpacked storage
561
* pos3: position of first eigenvalue of in sigx and sigz
564
for (i=0, pos=0, pos2=0, pos3=0; i<L;
565
pos += blck_szs[i]*(blck_szs[i]+1)/2,
566
pos2 += blck_szs[i]*blck_szs[i],
567
pos3 += blck_szs[i], i++){
569
lngth = blck_szs[i]*(blck_szs[i]+1)/2;
571
/* compute ith block of dF = \sum \delta x_i F_i,
573
F2C(dgemv)("N", &lngth, &m, &dbl1, F+sz+pos, &sz, dx, &int1,
576
/* scale dF as R'*dF*R, store in temp + lngth */
577
cngrncb(2, blck_szs[i], temp, R+pos2, temp+lngth, temp+2*lngth);
579
/* add (q/gap)*Lambda^(1/2) - Lambda^(-1/2) */
580
for (k=0, pos4=lngth; k<blck_szs[i]; pos4+=blck_szs[i]-k, k++)
581
temp[pos4] -= rho*sigx[pos3+k] + 1.0/sigx[pos3+k];
583
/* replace dF in temp by L^{-1}*dF*L^{-T},
584
* (L: cholesky factor of F, stored in X)
585
* and compute eigenvalues of L^{-1}*dF*L^{-T} */
586
F2C(dspgst)(&int1, "L", blck_szs+i, temp, X+pos, &info2);
588
sprintf(str,"Error in dspst, info = %d.\n", info2);
590
*info = -18; return 1;
592
/* temp has to be of size max_n*(max_n+1)+3*max_n */
593
F2C(dspev)("N", "L", blck_szs+i, temp, sigx+pos3, NULL, &int1,
594
temp+2*lngth, &info2);
596
sprintf(str,"Error in dspev, info = %d.\n", info2);
598
*info = -18; return 1;
601
/* dZ := R*((q/gap)*Lambda^(1/2) - Lambda^(-1/2) + R'*dF*R)*R' */
602
cngrncb(1, blck_szs[i], temp+lngth, R+pos2, dZ+pos,
605
/* copy ith block of dZ to temp */
606
F2C(dcopy)(&lngth, dZ+pos, &int1, temp, &int1);
608
/* copy ith block of Z to temp + lngth */
609
F2C(dcopy)(&lngth, Z+pos, &int1, temp+lngth, &int1);
611
/* sigz: generalized eigenvalues of (dZ,Z)
612
* required size of temp: 3*max_n + max_n*(max_n+1) */
613
F2C(dspgv)(&int1, "N", "L", blck_szs+i, temp, temp+lngth, sigz+pos3,
614
NULL, &int1, temp+2*lngth, &info2);
616
sprintf(str,"Error in dspgv, info = %d.\n", info2);Scistring(str);
617
*info = -18; return 1;
624
* compute feasible rectangle for plane search
627
maxpossigx = 0.0; minnegsigx = 0.0;
628
maxpossigz = 0.0; minnegsigz = 0.0;
629
for (i=0; i<n; i++) {
630
if ( sigx[i] > maxpossigx )
631
maxpossigx = sigx[i]; /* max pos eigenvalue in sigx */
632
else if ( sigx[i] < minnegsigx )
633
minnegsigx = sigx[i]; /* min neg eigenvalue in sigx */
634
if ( sigz[i] > maxpossigz )
635
maxpossigz = sigz[i]; /* max pos eigenvalue in sigz */
636
else if ( sigz[i] < minnegsigz )
637
minnegsigz = sigz[i]; /* min neg eigenvalue in sigz */
639
nrmx = F2C(dnrm2)(&n, sigx, &int1); /* norm of scaled dx */
640
nrmz = F2C(dnrm2)(&n, sigz, &int1); /* norm of scaled dZ */
641
nrmmax = MAX( nrmx, nrmz);
643
XdZ = inprd(F,dZ,L,blck_szs); /* Tr F0*dZ */
644
ZdX = F2C(ddot)(&m, c, &int1, dx, &int1); /* c^T*dx */
648
* check corners of feasible rectangle
651
dbl_epsilon = F2C(dlamch)("e");
652
if (nrmx > SIGTOL*nrmmax)
654
alphax = (minnegsigx < -dbl_epsilon) ? -1.0/minnegsigx : 0.0;
656
alphax = (maxpossigx > dbl_epsilon) ? -1.0/maxpossigx : 0.0;
659
if (nrmz > SIGTOL*nrmmax)
661
alphaz = (minnegsigz < -dbl_epsilon) ? -1.0/minnegsigz : 0.0;
663
alphaz = (maxpossigz > dbl_epsilon) ? -1.0/maxpossigz : 0.0;
666
newgap = gap + alphax*ZdX + alphaz*XdZ;
667
newu = ul[0] + alphax*ZdX;
668
newl = ul[1] - alphaz*XdZ;
670
if (newgap <= MAX(abstol, MINABSTOL)) *info = 2;
671
else if ( (newl > 0.0 && newgap <= reltol*newl) ||
672
(newu < 0.0 && newgap <= -reltol*newu) ) *info = 3;
673
else if ( reltol < 0.0 && newu <= tv ) *info = 4;
674
else if ( reltol < 0.0 && newl >= tv ) *info = 5;
675
else if ( *iters == maxiters ) *info = 1;
679
F2C(daxpy)(&m, &alphax, dx, &int1, x, &int1);
680
F2C(daxpy)(&sz, &alphaz, dZ, &int1, Z, &int1);
681
gap = newgap; ul[0] = newu; ul[1] = newl;
682
sprintf(str,"% 13.2e % 12.2e %10.2e \n", ul[0], ul[1], gap);
691
* minimize phi(alphax,alphaz) =
692
* q*log(dual_gap + alphax*c^T*dx + alphaz* Tr F_0 dZ)
693
* - sum log (1+alphax*sigx_i) - sum log (1+alphaz*sigz)
696
alphax = 0.0; alphaz = 0.0; lambda_ls = 1.0;
698
if (nrmx > SIGTOL*nrmmax)
699
if (nrmz > SIGTOL*nrmmax) /* compute primal and dual steps */
700
while ( lambda_ls > 1e-4 ) {
702
/* compute 1st and 2nd derivatives of phi */
703
rho = q/(gap + alphax*ZdX + alphaz*XdZ);
704
gradx = rho*ZdX; hessx = 0.0;
705
gradz = rho*XdZ; hessz = 0.0;
707
gradx -= sigx[i] / (1.0+alphax*sigx[i]);
708
hessx += SQR( sigx[i] / (1.0+alphax*sigx[i]) );
709
gradz -= sigz[i] / (1.0+alphaz*sigz[i]);
710
hessz += SQR( sigz[i] / (1.0+alphaz*sigz[i]) );
714
dalphax = -gradx/hessx; dalphaz = -gradz/hessz;
715
lambda_ls = sqrt( SQR(gradx)/hessx + SQR(gradz)/hessz );
716
alphax += (lambda_ls > 0.25) ?
717
dalphax/(1.0+lambda_ls) : dalphax;
718
alphaz += (lambda_ls > 0.25) ?
719
dalphaz/(1.0+lambda_ls) : dalphaz;
723
else while ( lambda_ls > 1e-4 ) { /* primal step only */
725
/* compute 1st and 2nd derivatives of phi */
726
rho = q/(gap + alphax*ZdX);
727
gradx = rho*ZdX; hessx = 0.0;
729
gradx -= sigx[i] / (1.0+alphax*sigx[i]);
730
hessx += SQR( sigx[i] / (1.0+alphax*sigx[i]) );
734
dalphax = -gradx/hessx;
735
lambda_ls = fabs(gradx)/sqrt(hessx);
736
alphax += (lambda_ls > 0.25) ?
737
dalphax/(1.0+lambda_ls) : dalphax;
741
else if (nrmz > SIGTOL*nrmmax) /* dual step only */
742
while ( lambda_ls > 1e-4 ) {
744
/* compute 1st and 2nd derivatives of phi */
745
rho = q/(gap + alphaz*XdZ);
746
gradz = rho*XdZ; hessz = 0.0;
748
gradz -= sigz[i] / (1.0+alphaz*sigz[i]);
749
hessz += SQR( sigz[i] / (1.0+alphaz*sigz[i]) );
753
dalphaz = -gradz/hessz;
754
lambda_ls = fabs(gradz)/sqrt(hessz);
755
alphaz += (lambda_ls > 0.25) ?
756
dalphaz/(1.0+lambda_ls) : dalphaz;
762
F2C(daxpy)(&m, &alphax, dx, &int1, x, &int1);
763
F2C(daxpy)(&sz, &alphaz, dZ, &int1, Z, &int1);
767
return -1; /* should never happen */
added
removed
routines/optim/sp.c
