parallel_loqo.c

支持向量分类算法在linux操作系统下的是实现

      /* 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;
}