gamma_z   = malloc(n*sizeof(double));  gamma_s   = malloc(n*sizeof(double));  hat_nu    = malloc(n*sizeof(double));  hat_tau   = malloc(n*sizeof(double));  delta_x   = malloc(n*sizeof(double));  delta_y   = malloc(m*sizeof(double));  delta_s   = malloc(n*sizeof(double));  delta_z   = malloc(n*sizeof(double));  delta_g   = malloc(n*sizeof(double));  delta_t   = malloc(n*sizeof(double));  d         = malloc(n*sizeof(double));  /* pointers into the external variables */  x = primal;			/* n */  g = x + n;			/* n */  t = g + n;			/* n */  y = dual;			/* m */  z = y + m;			/* n */  s = z + n;			/* n */  /* initial settings */  b_plus_1 = 1;  c_plus_1 = 0;  for (i=0; i<n; i++) c_plus_1 += c[i];  /* get diagonal terms */  for (i=0; i<n; i++) diag_h_x[i] = h_x[(n+1)*i];     /* starting point */  if (restart == 1) {				/* x, y already preset */    for (i=0; i<n; i++) {	/* compute g, t for primal feasibility */      g[i] = max(ABS(x[i] - l[i]), bound);      t[i] = max(ABS(u[i] - x[i]), bound);     }    matrix_vector(n, h_x, x, h_dot_x); /* h_dot_x = h_x * x */    for (i=0; i<n; i++) {	/* sigma is a dummy variable to calculate z, s */      sigma[i] = c[i] + h_dot_x[i];      for (j=0; j<m; j++)	sigma[i] -= a[n*j + i] * y[j];      if (sigma[i] > 0) {	s[i] = bound;	z[i] = sigma[i] + bound;      }      else {	s[i] = bound - sigma[i];	z[i] = bound;      }    }  }  else {			/* use default start settings */    for (i=0; i<m; i++)      for (j=i; j<m; j++)	h_y[i*m + j] = (i==j) ? 1 : 0;        for (i=0; i<n; i++) {      c_x[i] = c[i];      h_x[(n+1)*i] += 1;    }        for (i=0; i<m; i++)      c_y[i] = b[i];        /* and solve the system [-H_x A'; A H_y] [x, y] = [c_x; c_y] */    solve_reduced(n, m, h_x, h_y, a, x, y, c_x, c_y, workspace,		  PREDICTOR);        /* initialize the other variables */    for (i=0; i<n; i++) {      g[i] = max(ABS(x[i] - l[i]), bound);      z[i] = max(ABS(x[i]), bound);      t[i] = max(ABS(u[i] - x[i]), bound);       s[i] = max(ABS(x[i]), bound);     }  }  for (i=0, mu=0; i<n; i++)    mu += z[i] * g[i] + s[i] * t[i];  mu = mu / (2*n);  /* the main loop */  if (verb >= STATUS) {    printf("counter | pri_inf  | dual_inf  | pri_obj   | dual_obj  | ");    printf("sigfig | alpha  | nu \n");    printf("-------------------------------------------------------");    printf("---------------------------\n");  }    while (status == STILL_RUNNING) {    /* predictor */        /* put back original diagonal values */    for (i=0; i<n; i++)       h_x[(n+1) * i] = diag_h_x[i];    matrix_vector(n, h_x, x, h_dot_x); /* compute h_dot_x = h_x * x */    for (i=0; i<m; i++) {      rho[i] = b[i];      for (j=0; j<n; j++)	rho[i] -= a[n*i + j] * x[j];    }        for (i=0; i<n; i++) {      nu[i] = l[i] - x[i] + g[i];      tau[i] = u[i] - x[i] - t[i];      sigma[i] = c[i] - z[i] + s[i] + h_dot_x[i];      for (j=0; j<m; j++)	sigma[i] -= a[n*j + i] * y[j];      gamma_z[i] = - z[i];      gamma_s[i] = - s[i];    }    /* instrumentation */    x_h_x = 0;    primal_inf = 0;    dual_inf = 0;        for (i=0; i<n; i++) {      x_h_x += h_dot_x[i] * x[i];      primal_inf += sqr(tau[i]);      primal_inf += sqr(nu[i]);      dual_inf += sqr(sigma[i]);    }    for (i=0; i<m; i++)       primal_inf += sqr(rho[i]);    primal_inf = sqrt(primal_inf)/b_plus_1;    dual_inf = sqrt(dual_inf)/c_plus_1;    primal_obj = 0.5 * x_h_x;    dual_obj = -0.5 * x_h_x;    for (i=0; i<n; i++) {      primal_obj += c[i] * x[i];      dual_obj += l[i] * z[i] - u[i] * s[i];    }    for (i=0; i<m; i++)      dual_obj += b[i] * y[i];    sigfig = log10(ABS(primal_obj) + 1) -             log10(ABS(primal_obj - dual_obj));    sigfig = max(sigfig, 0);		       /* the diagnostics - after we computed our results we will       analyze them */    if (counter > counter_max) status = ITERATION_LIMIT;    if (sigfig  > sigfig_max)  status = OPTIMAL_SOLUTION;    if (primal_inf > 10e100)   status = PRIMAL_INFEASIBLE;    if (dual_inf > 10e100)     status = DUAL_INFEASIBLE;    if ((primal_inf > 10e100) & (dual_inf > 10e100)) status = PRIMAL_AND_DUAL_INFEASIBLE;    if (ABS(primal_obj) > 10e100) status = PRIMAL_UNBOUNDED;    if (ABS(dual_obj) > 10e100) status = DUAL_UNBOUNDED;    /* write some nice routine to enforce the time limit if you       _really_ want, however it's quite useless as you can compute       the time from the maximum number of iterations as every       iteration costs one cholesky decomposition plus a couple of        backsubstitutions */    /* generate report */    if ((verb >= FLOOD) | ((verb == STATUS) & (status != 0)))      printf("%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n",	     counter, primal_inf, dual_inf, primal_obj, dual_obj,	     sigfig, alfa, mu);    counter++;    if (status == 0) {		/* we may keep on going, otherwise				   it'll cost one loop extra plus a				   messed up main diagonal of h_x */      /* intermediate variables (the ones with hat) */      for (i=0; i<n; i++) {	hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];	hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];	/* diagonal terms */	d[i] = z[i] / g[i] + s[i] / t[i];      }            /* initialization before the cholesky solver */      for (i=0; i<n; i++) {	h_x[(n+1)*i] = diag_h_x[i] + d[i];	c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - 	  s[i] * hat_tau[i] / t[i];      }      for (i=0; i<m; i++) {	c_y[i] = rho[i];	for (j=i; j<m; j++) 	  h_y[m*i + j] = 0;      }            /* and do it */      solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,		    PREDICTOR);            for (i=0; i<n; i++) {	/* backsubstitution */	delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];	delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];		delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];	delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];		/* central path (corrector) */	gamma_z[i] = mu / g[i] - z[i] - delta_z[i] * delta_g[i] / g[i];	gamma_s[i] = mu / t[i] - s[i] - delta_s[i] * delta_t[i] / t[i];		/* (some more intermediate variables) the hat variables */	hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];	hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];		/* initialization before the cholesky */	c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - s[i] * hat_tau[i] / t[i];      }            for (i=0; i<m; i++) {	/* comput c_y and rho */	c_y[i] = rho[i];	for (j=i; j<m; j++)	  h_y[m*i + j] = 0;      }            /* and do it */      solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,		    CORRECTOR);            for (i=0; i<n; i++) {	/* backsubstitution */	delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];	delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];		delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];	delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];      }            alfa = -1;      for (i=0; i<n; i++) {	alfa = min(alfa, delta_g[i]/g[i]);	alfa = min(alfa, delta_t[i]/t[i]);	alfa = min(alfa, delta_s[i]/s[i]);	alfa = min(alfa, delta_z[i]/z[i]);      }      alfa = (margin - 1) / alfa;            /* compute mu */      for (i=0, mu=0; i<n; i++)	mu += z[i] * g[i] + s[i] * t[i];      mu = mu / (2*n);      mu = mu * sqr((alfa - 1) / (alfa + 10));            for (i=0; i<n; i++) {	x[i] += alfa * delta_x[i];	g[i] += alfa * delta_g[i];	t[i] += alfa * delta_t[i];	z[i] += alfa * delta_z[i];	s[i] += alfa * delta_s[i];      }            for (i=0; i<m; i++) 	y[i] += alfa * delta_y[i];    }  }  if ((status == 1) && (verb >= STATUS)) {    printf("----------------------------------------------------------------------------------\n");    printf("optimization converged\n");  }    /* free memory */  free(workspace);  free(diag_h_x);  free(h_y);  free(c_x);  free(c_y);  free(h_dot_x);    free(rho);  free(nu);  free(tau);  free(sigma);  free(gamma_z);  free(gamma_s);    free(hat_nu);  free(hat_tau);      free(delta_x);  free(delta_y);  free(delta_s);  free(delta_z);  free(delta_g);  free(delta_t);      free(d);  /* and return to sender */  return status;}