29 #include <atlas/clapack.h>
47 SG_SDEBUG(
"terminating optimizer - %s\n", error_text);
68 for (int32_t i=0; i<n; i++)
70 for (int32_t j=0; j<n; j++)
75 int result=clapack_dpotrf(CblasRowMajor, CblasUpper, (
int) n, a2, (
int) n);
77 for (int32_t i=0; i<n; i++)
80 for (int32_t i=0; i<n; i++)
82 for (int32_t j=i+1; j<n; j++)
89 SG_SDEBUG(
"Choldc failed, matrix not positive definite\n");
98 void nrerror(
char error_text[]);
102 for (i = 0; i < n; i++)
104 for (j = i; j < n; j++)
108 for (k=i-1; k>=0; k--)
109 sum -= a[n*i + k]*a[n*j + k];
115 SG_SDEBUG(
"Choldc failed, matrix not positive definite");
124 a[n*j + i] = sum/p[i];
138 for (i=0; i<n; i++) {
140 for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
144 for (i=n-1; i>=0; i--) {
146 for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
162 for (i=0; i<n; i++) {
164 for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
175 for (i=n-1; i>=0; i--) {
177 for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
236 h_y[m*i + j] += t_a[n*j + k] * t_a[n*i + k];
244 for (i=0; i<m; i++) {
247 t_y[i] += t_a[i*n + j] * t_c[j];
250 cholsb(h_y, m, p_y, t_y, x_y);
252 for (i=0; i<n; i++) {
255 t_c[i] += t_a[j*n + i] * x_y[j];
274 for (i=0; i<n; i++) {
275 y[i] = m[(n+1) * i] * x[i];
278 y[i] += m[i + n*j] * x[j];
280 for (j=i+1; j<n; j++)
281 y[i] += m[n*i + j] * x[j];
403 for (i=0; i<n; i++) c_plus_1 += c[i];
406 for (i=0; i<n; i++) diag_h_x[i] = h_x[(n+1)*i];
411 for (i=0; i<n; i++) {
418 for (i=0; i<n; i++) {
419 sigma[i] = c[i] + h_dot_x[i];
421 sigma[i] -= a[n*j + i] * y[j];
425 z[i] = sigma[i] + bound;
428 s[i] = bound - sigma[i];
436 h_y[i*m + j] = (i==j) ? 1 : 0;
438 for (i=0; i<n; i++) {
447 solve_reduced(n, m, h_x, h_y, a, x, y, c_x, c_y, workspace,
451 for (i=0; i<n; i++) {
459 for (i=0, mu=0; i<n; i++)
460 mu += z[i] * g[i] + s[i] * t[i];
465 SG_SDEBUG(
"counter | pri_inf | dual_inf | pri_obj | dual_obj | ");
467 SG_SDEBUG(
"-------------------------------------------------------");
468 SG_SDEBUG(
"---------------------------\n");
476 h_x[(n+1) * i] = diag_h_x[i];
480 for (i=0; i<m; i++) {
483 rho[i] -= a[n*i + j] * x[j];
486 for (i=0; i<n; i++) {
487 nu[i] = l[i] - x[i] + g[i];
488 tau[i] = u[i] - x[i] - t[i];
490 sigma[i] = c[i] - z[i] + s[i] + h_dot_x[i];
492 sigma[i] -= a[n*j + i] * y[j];
503 for (i=0; i<n; i++) {
504 x_h_x += h_dot_x[i] * x[i];
511 primal_inf = sqrt(primal_inf)/b_plus_1;
512 dual_inf = sqrt(dual_inf)/c_plus_1;
514 primal_obj = 0.5 * x_h_x;
515 dual_obj = -0.5 * x_h_x;
516 for (i=0; i<n; i++) {
517 primal_obj += c[i] * x[i];
518 dual_obj += l[i] * z[i] - u[i] * s[i];
521 dual_obj += b[i] * y[i];
545 if ((verb >=
FLOOD) | ((verb ==
STATUS) & (status != 0)))
546 SG_SDEBUG(
"%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n",
547 counter, primal_inf, dual_inf, primal_obj, dual_obj,
556 for (i=0; i<n; i++) {
557 hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
558 hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
560 d[i] = z[i] / g[i] + s[i] / t[i];
564 for (i=0; i<n; i++) {
565 h_x[(n+1)*i] = diag_h_x[i] + d[i];
566 c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] -
567 s[i] * hat_tau[i] / t[i];
569 for (i=0; i<m; i++) {
576 if (!
solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,
PREDICTOR))
582 for (i=0; i<n; i++) {
584 delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
585 delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
587 delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
588 delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
591 gamma_z[i] = mu / g[i] - z[i] - delta_z[i] * delta_g[i] / g[i];
592 gamma_s[i] = mu / t[i] - s[i] - delta_s[i] * delta_t[i] / t[i];
595 hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
596 hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
599 c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - s[i] * hat_tau[i] / t[i];
602 for (i=0; i<m; i++) {
609 solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,
612 for (i=0; i<n; i++) {
614 delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
615 delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
617 delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
618 delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
622 for (i=0; i<n; i++) {
628 alfa = (margin - 1) / alfa;
631 for (i=0, mu=0; i<n; i++)
632 mu += z[i] * g[i] + s[i] * t[i];
634 mu = mu *
CMath::sq((alfa - 1) / (alfa + 10));
636 for (i=0; i<n; i++) {
637 x[i] += alfa * delta_x[i];
638 g[i] += alfa * delta_g[i];
639 t[i] += alfa * delta_t[i];
640 z[i] += alfa * delta_z[i];
641 s[i] += alfa * delta_s[i];
645 y[i] += alfa * delta_y[i];
650 if ((status == 1) && (verb >=
STATUS)) {
651 SG_SDEBUG(
"----------------------------------------------------------------------------------\n");