/* MPI Barrier to wait for changes in x. */      info = MPI_Barrier(comm); CheckError(info);      /* Initialize the other variables besides x and y. */      for(i=0; i<local_n; ++i)	{	  idx = i*local_stride;	  loc.g[idx] = max(ABS(loc.x[idx] - loc.l[idx]), bound);	  loc.z[idx] = max(ABS(loc.x[idx]), bound);	  loc.t[idx] = max(ABS(loc.u[idx] - loc.x[idx]), bound);	  loc.s[idx] = max(ABS(loc.x[idx]), bound);	}    }  /* MPI Barrier to wait for local changes in s,t,z and g. */  info = MPI_Barrier(comm); CheckError(info);  /* Set mu <- (s^T * t + z^T * g) / (2*n) */  info = PLA_Dot(*z, *g, mu); CheckError(info);  info = PLA_Dot(*s, *t, tmp); CheckError(info);  info = MPI_Barrier(comm);  *loc.mu = (*loc.mu + (*loc.tmp)) / (2*n);/*   printf("Setting initial values done.\n"); fflush(stdout); */  /* The main loop. */  if(verb >= STATUS && rank == 0)    {      printf("counter | pri_inf  | dual_inf  | pri_obj   | dual_obj  | ");      printf("sigfig | alpha  | nu \n");      printf("-------------------------------------------------------");      printf("---------------------------\n");    }  while(status == STILL_RUNNING)    {      /* Predictor step. */      /* Put back original diagonal values. */      PLA_Set_diagonal(&h_x, diag_h_x);      /* Set h_dot_x <- h_x * x *//*       PLA_Copy_vector(*x, h_dot_x); */      info = PLA_Obj_set_to_zero(h_dot_x); CheckError(info);/*        info = PLA_Symv(PLA_UPPER_TRIANGULAR, plus_one, h_x, *x, plus_one, h_dot_x);  */      info = PLA_Gemv(PLA_NO_TRANS, plus_one, h_x_copy, *x, plus_one, 		      h_dot_x); CheckError(info);      /* Set rho <- b - A * x */      PLA_Copy_vector(b, rho);      info = PLA_Gemv(PLA_NO_TRANS, minus_one, a, *x, plus_one, rho);       CheckError(info);      /* Set nu <- l - x + g */      PLA_Copy_vector(l, nu);      info = PLA_Axpy(minus_one, *x, nu); CheckError(info);      info = PLA_Axpy(plus_one, *g, nu); CheckError(info);      /* Set tau <- u - x - t */      PLA_Copy_vector(u, tau);      info = PLA_Axpy(minus_one, *x, tau); CheckError(info);      info = PLA_Axpy(minus_one, *t, tau); CheckError(info);      /* Set sigma <- c - z + s + h_dot_x - A^T * y */      PLA_Copy_vector(c, sigma);      info = PLA_Axpy(minus_one, *z, sigma); CheckError(info);      info = PLA_Axpy(plus_one, *s, sigma); CheckError(info);      info = PLA_Axpy(plus_one, h_dot_x, sigma); CheckError(info);      info = PLA_Gemv(PLA_TRANS, minus_one, a, *y, plus_one, sigma);      /* Set gamma_z <- -z, gamma_s <- -s. */      info = PLA_Obj_set_to_zero(gamma_s); CheckError(info);      info = PLA_Axpy(minus_one, *s, gamma_s); CheckError(info);      info = PLA_Obj_set_to_zero(gamma_z); CheckError(info);      info = PLA_Axpy(minus_one, *z, gamma_z); CheckError(info);/*       print_vector(h_dot_x); *//*       print_vector(rho); *//*       print_vector(nu); *//*       print_vector(tau); *//*       print_vector(sigma); */      /* Instrumentation */      info = PLA_Dot(h_dot_x, *x, x_h_x); CheckError(info);      {	/* Compute 	   primal_inf <- sqrt(tau^T * tau + nu^T * nu + rho^T * rho) 	   / b_plus_1. */	info = PLA_Dot(tau, tau, primal_inf); CheckError(info);	info = PLA_Dot(nu, nu, tmp); CheckError(info);	info = MPI_Barrier(comm);	*loc.primal_inf += *loc.tmp;	info = PLA_Dot(rho, rho, tmp); CheckError(info);	info = MPI_Barrier(comm);	*loc.primal_inf += *loc.tmp;	*loc.primal_inf = sqrt(*loc.primal_inf)/(*loc.b_plus_1);      }      {	/* Compute 	   dual_inf <- sqrt(sigma^T * sigma) / c_plus_1. */	info = PLA_Dot(sigma, sigma, dual_inf); CheckError(info);	info = MPI_Barrier(comm);	*loc.dual_inf = sqrt(*loc.dual_inf)/(*loc.c_plus_1);      }      {	/* Compute	   primal_obj <- 0.5 * x_h_x + c^T * x. */	info = PLA_Dot(c, *x, primal_obj); CheckError(info);	info = MPI_Barrier(comm);	*loc.primal_obj += 0.5*(*loc.x_h_x);      }      {	/* Compute	   dual_obj <- -0.5 * x_h_x + l^T * z - u^T * s + b^T * y. */	info = PLA_Dot(l, *z, dual_obj); CheckError(info);	info = PLA_Dot(b, *y, tmp); CheckError(info);	info = MPI_Barrier(comm);	*loc.dual_obj += *loc.tmp;	info = PLA_Dot(u, *s, tmp); CheckError(info);	info = MPI_Barrier(comm);	*loc.dual_obj -= *loc.tmp;	*loc.dual_obj -= 0.5*(*loc.x_h_x);      }      sigfig = log10(ABS(*loc.primal_obj) + 1) -	log10(ABS(*loc.primal_obj - (*loc.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 (*loc.primal_inf > 10e100)   status = PRIMAL_INFEASIBLE;      if (*loc.dual_inf > 10e100)     status = DUAL_INFEASIBLE;      if (*loc.primal_inf > 10e100 && *loc.dual_inf > 10e100) 	status = PRIMAL_AND_DUAL_INFEASIBLE;      if (ABS(*loc.primal_obj) > 10e100) status = PRIMAL_UNBOUNDED;      if (ABS(*loc.dual_obj) > 10e100) status = DUAL_UNBOUNDED;            /* Generate report */      if ((verb >= FLOOD) | ((verb == STATUS) & (status != 0)) && rank == 0)	printf("%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n",	       counter, *loc.primal_inf, *loc.dual_inf, *loc.primal_obj, 	       *loc.dual_obj, sigfig, *loc.alfa, *loc.mu);      counter++;      if(status == 0){	/* Set intermediary hat variables. */	for(i=0; i<local_n; ++i)	  {	    idx = i*local_stride;	    loc.hat_nu[idx] = loc.nu[idx] + loc.g[idx] * 	      loc.gamma_z[idx] / loc.z[idx];	    loc.hat_tau[idx] = loc.tau[idx] - loc.t[idx] *	      loc.gamma_s[idx] / loc.s[idx];	    loc.d[idx] = loc.z[idx] / loc.g[idx] + loc.s[idx] / loc.t[idx];	    loc.c_x[idx] = loc.sigma[idx] - loc.z[idx] * loc.hat_nu[idx] /	      loc.g[idx] - loc.s[idx] * loc.hat_tau[idx] / loc.t[idx];	  }	info = MPI_Barrier(comm); CheckError(info);	/* Set h_x <- h_x + diag(z^-1 * g + s * t^-1). */	info = PLA_Copy(h_x_copy, h_x); CheckError(info);	info = PLA_Axpy(plus_one, diag_h_x, d); CheckError(info);	PLA_Set_diagonal(&h_x, d);	/* Set c_y <- rho. */	PLA_Copy_vector(rho, c_y);	/* Set h_y <-  0. */	info = PLA_Obj_set_to_zero(h_y); CheckError(info);		/* Compute predictor step */	parallel_solve_reduced(h_x, h_y, a, c_x, c_y, &delta_x, &delta_y, &y1t, 			       PREDICTOR);	info = MPI_Barrier(comm); CheckError(info);	/* Do backsubstitution */	for(i=0; i<local_n; ++i)	  {	    idx = i*local_stride;	    loc.delta_s[idx] = loc.s[idx] * 	      (loc.delta_x[idx] - loc.hat_tau[idx]) / loc.t[idx];	    loc.delta_z[idx] = loc.z[idx] *	      (loc.hat_nu[idx] - loc.delta_x[idx]) / loc.g[idx];	    loc.delta_g[idx] = loc.g[idx] *	      (loc.gamma_z[idx] - loc.delta_z[idx]) / loc.z[idx];	    loc.delta_t[idx] = loc.t[idx] *	      (loc.gamma_s[idx] - loc.delta_s[idx]) / loc.s[idx];	    	    /* Central path (corrector) */	    loc.gamma_z[idx] = *loc.mu / loc.g[idx] - loc.z[idx] -	      loc.delta_z[idx] * loc.delta_g[idx] / loc.g[idx];	    loc.gamma_s[idx] = *loc.mu / loc.t[idx] - loc.s[idx] -	      loc.delta_s[idx] * loc.delta_t[idx] / loc.t[idx];	    	    /* The hat variables. */	    loc.hat_nu[idx] = loc.nu[idx] + loc.g[idx] * 	      loc.gamma_z[idx] / loc.z[idx];	    loc.hat_tau[idx] = loc.tau[idx] - loc.t[idx] *	      loc.gamma_s[idx] / loc.s[idx];	    	    loc.c_x[idx] = loc.sigma[idx] - loc.z[idx] * loc.hat_nu[idx] /	      loc.g[idx] - loc.s[idx] * loc.hat_tau[idx] / loc.t[idx];	  }	info = MPI_Barrier(comm); CheckError(info);	/* Set c_y */	PLA_Copy_vector(rho, c_y);	/* Compute corrector step */	parallel_solve_reduced(h_x, h_y, a, c_x, c_y, &delta_x, &delta_y, &y1t, 			       CORRECTOR);		info = MPI_Barrier(comm); CheckError(info);	/* Backsubstitution. */	for(i=0; i<local_n; ++i)	  {	    idx = i*local_stride;	    loc.delta_s[idx] = loc.s[idx] * 	      (loc.delta_x[idx] - loc.hat_tau[idx]) / loc.t[idx];	    loc.delta_z[idx] = loc.z[idx] *	      (loc.hat_nu[idx] - loc.delta_x[idx]) / loc.g[idx];	    loc.delta_g[idx] = loc.g[idx] *	      (loc.gamma_z[idx] - loc.delta_z[idx]) / loc.z[idx];	    loc.delta_t[idx] = loc.t[idx] *	      (loc.gamma_s[idx] - loc.delta_s[idx]) / loc.s[idx];	  }	/* Update alfa. */	for(i=0; i<local_n; ++i)	  {	    idx = i*local_stride;	    loc.max_ratio[idx] = -loc.delta_g[idx] / loc.g[idx];	    next_ratio = -loc.delta_t[idx] / loc.t[idx];	    if(loc.max_ratio[idx] < next_ratio)	      loc.max_ratio[idx] = next_ratio;	    next_ratio = -loc.delta_s[idx] / loc.s[idx];	    if(loc.max_ratio[idx] < next_ratio)	      loc.max_ratio[idx] = next_ratio;	    next_ratio = -loc.delta_z[idx] / loc.z[idx];	    if(loc.max_ratio[idx] < next_ratio)	      loc.max_ratio[idx] = next_ratio;	  }	info = MPI_Barrier(comm); CheckError(info);	/* Note: There is a bug in the plapack documentation!	 * Correct calling sequence is:	 * int PLA_Iamax( PLA_Obj x, PLA_Obj k, PLA_Obj xmax). */	info = PLA_Iamax(max_ratio, max_idx, alfa); CheckError(info);	info = MPI_Barrier(comm); CheckError(info);	*loc.alfa = -max(*loc.alfa, 1);	*loc.alfa = (margin - 1) / (*loc.alfa);		/* Set mu <- (s^T * t + z^T * g) / (2*n) * (	   (alfa-1) / (alfa + 10))^2. */	info = PLA_Dot(*z, *g, mu); CheckError(info);	info = PLA_Dot(*s, *t, tmp); CheckError(info);	info = MPI_Barrier(comm); CheckError(info);	*loc.mu = (*loc.mu + (*loc.tmp)) / (2*n);	*loc.mu = *loc.mu * sqr((*loc.alfa - 1) / (*loc.alfa + 10));	/* Update variables. */	info = PLA_Axpy(alfa, delta_x, *x); CheckError(info);	info = PLA_Axpy(alfa, delta_g, *g); CheckError(info);	info = PLA_Axpy(alfa, delta_t, *t); CheckError(info);	info = PLA_Axpy(alfa, delta_z, *z); CheckError(info);	info = PLA_Axpy(alfa, delta_s, *s); CheckError(info);	info = PLA_Axpy(alfa, delta_y, *y); CheckError(info);      } /* if(status == 0) */    } /* while(status == STILL_RUNNING) */  if(rank == 0 && status == 1 && verb >= STATUS)    {      printf("-----------------------------------------");      printf("-----------------------------------------\n");      printf("optimization converged\n");    }  info = MPI_Barrier(comm); CheckError(info);  /* Compute dist <- c + h_dot_x - A^T*y */  PLA_Copy_vector(h_dot_x, *dist);  info = PLA_Gemv(PLA_TRANS, minus_one, a, *y, plus_one, *dist);  CheckError(info);  info = PLA_Axpy(plus_one, c, *dist); CheckError(info);  info = MPI_Barrier(comm); CheckError(info);  /* Deallocation */  info = PLA_Obj_free(&diag_h_x); CheckError(info);    info = PLA_Obj_free(&h_y); CheckError(info);    info = PLA_Obj_free(&y1t); CheckError(info);    info = PLA_Obj_free(&c_x); CheckError(info);  info = PLA_Obj_free(&c_y); CheckError(info);  info = PLA_Obj_free(&h_dot_x); CheckError(info);  info = PLA_Obj_free(&rho); CheckError(info);  info = PLA_Obj_free(&nu); CheckError(info);  info = PLA_Obj_free(&tau); CheckError(info);  info = PLA_Obj_free(&sigma); CheckError(info);  info = PLA_Obj_free(&gamma_z); CheckError(info);  info = PLA_Obj_free(&gamma_s); CheckError(info);  info = PLA_Obj_free(&hat_nu); CheckError(info);  info = PLA_Obj_free(&hat_tau); CheckError(info);  info = PLA_Obj_free(&delta_x); CheckError(info);  info = PLA_Obj_free(&delta_y); CheckError(info);  info = PLA_Obj_free(&delta_s); CheckError(info);  info = PLA_Obj_free(&delta_z); CheckError(info);  info = PLA_Obj_free(&delta_g); CheckError(info);  info = PLA_Obj_free(&delta_t); CheckError(info);  info = PLA_Obj_free(&d); CheckError(info);  info = PLA_Obj_free(&ones_n); CheckError(info);  info = PLA_Obj_free(&ones_m); CheckError(info);  info = PLA_Obj_free(&h_x_copy); CheckError(info);  info = PLA_Obj_free(&plus_one); CheckError(info);  info = PLA_Obj_free(&minus_one); CheckError(info);  info = PLA_Obj_free(&b_plus_1); CheckError(info);  info = PLA_Obj_free(&c_plus_1); CheckError(info);  info = PLA_Obj_free(&x_h_x); CheckError(info);  info = PLA_Obj_free(&primal_inf); CheckError(info);  info = PLA_Obj_free(&dual_inf); CheckError(info);  info = PLA_Obj_free(&primal_obj); CheckError(info);  info = PLA_Obj_free(&dual_obj); CheckError(info);  info = PLA_Obj_free(&mu); CheckError(info);  info = PLA_Obj_free(&alfa); CheckError(info);  info = PLA_Obj_free(&tmp); CheckError(info);  return status;}