parallel_loqo.c

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

}

/* Struct for local buffer access. */
struct local_buffer_address{
  LOQOfloat *c;
  LOQOfloat *b;

  LOQOfloat *x;
  LOQOfloat *l;
  LOQOfloat *u;
  LOQOfloat *g;
  LOQOfloat *t;
  LOQOfloat *s;
  LOQOfloat *z;

  LOQOfloat *delta_x;
  LOQOfloat *delta_g;
  LOQOfloat *delta_t;
  LOQOfloat *delta_z;
  LOQOfloat *delta_s;

  LOQOfloat *sigma;
  LOQOfloat *nu;
  LOQOfloat *tau;
  LOQOfloat *hat_nu;
  LOQOfloat *hat_tau;
  LOQOfloat *gamma_z;
  LOQOfloat *gamma_s;

  LOQOfloat *c_plus_1;
  LOQOfloat *b_plus_1;
  LOQOfloat *primal_inf;
  LOQOfloat *dual_inf;
  LOQOfloat *primal_obj;
  LOQOfloat *dual_obj;
  LOQOfloat *x_h_x;
  LOQOfloat *tmp;
  LOQOfloat *mu;
  LOQOfloat *alfa;
  LOQOfloat *max_ratio;

  LOQOfloat *d;
  LOQOfloat *c_x;
  LOQOfloat *c_y;
};
   
int parallel_loqo(PLA_Obj c, PLA_Obj h_x, PLA_Obj a, PLA_Obj b, PLA_Obj l,
		  PLA_Obj u, PLA_Obj *x, PLA_Obj *g, PLA_Obj *t, PLA_Obj *y,
		  PLA_Obj *z, PLA_Obj *s, PLA_Obj *dist, int verb, LOQOfloat sigfig_max, 
		  int counter_max, LOQOfloat margin, LOQOfloat bound, int restart)
{
  /* To be allocated. */
  PLA_Obj diag_h_x = NULL, h_y = NULL, y1t = NULL, c_x = NULL, c_y = NULL; 
  PLA_Obj h_dot_x = NULL, rho = NULL, nu = NULL, tau = NULL, sigma = NULL;
  PLA_Obj gamma_z = NULL, gamma_s = NULL, hat_nu = NULL, hat_tau = NULL;
  PLA_Obj delta_x = NULL, delta_y = NULL, delta_s = NULL, delta_z = NULL;
  PLA_Obj delta_g = NULL, delta_t = NULL, d = NULL, ones_n = NULL;
  PLA_Obj ones_m = NULL, max_ratio = NULL;
  PLA_Obj h_x_copy = NULL;
  
  /* Multiscalars */
  PLA_Obj b_plus_1 = NULL, c_plus_1 = NULL, x_h_x = NULL, primal_inf = NULL;
  PLA_Obj dual_inf = NULL, primal_obj = NULL, dual_obj = NULL, mu = NULL;
  PLA_Obj alfa = NULL, plus_one = NULL, minus_one = NULL, tmp = NULL;
  PLA_Obj max_idx = NULL;

  /* Local buffer stuff */
  void *addr_buf;
  struct local_buffer_address loc;
  int idx, local_n, local_stride, local_m, local_stride2;

  /* Instrumentation etc. */
  LOQOfloat sigfig = 0.0;
  LOQOfloat next_ratio = 0.0;
  int counter = 0;
  int status = STILL_RUNNING;
  int i,n,m,info;
  PLA_Template templ = NULL;
  MPI_Comm comm;
  int rank, numnodes;

  info = PLA_Set_error_checking(FALSE,FALSE,FALSE,FALSE);
  CheckError(info);

  /* Extract template and MPI information */
  info = PLA_Obj_template(h_x, &templ); CheckError(info);
  info = PLA_Temp_comm_all_info(templ, &comm, &rank, &numnodes); 
  CheckError(info);
  info = PLA_Obj_global_length(h_x, &n); CheckError(info);
  info = PLA_Obj_global_length(b, &m); CheckError(info);

/*   print_vector(c); */
/*   print_matrix(h_x); */
/*   print_matrix(a); */
/*   print_vector(b); */
/*   print_vector(l); */
/*   print_vector(u); */

/*   printf("Allocating objs..."); fflush(stdout); */
  /* Allocation */
  info = PLA_Mvector_create_conf_to(c, 1, &diag_h_x); CheckError(info);
  info = PLA_Matrix_create(MPIfloat, m, m, templ, PLA_ALIGN_FIRST, 
			   PLA_ALIGN_FIRST, &h_y); CheckError(info);
  info = PLA_Matrix_create(MPIfloat, m, n, templ, PLA_ALIGN_FIRST,
			   PLA_ALIGN_FIRST, &y1t); CheckError(info);
  info = PLA_Matrix_create(MPIfloat, n, n, templ, PLA_ALIGN_FIRST,
			   PLA_ALIGN_FIRST, &h_x_copy); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &c_x); CheckError(info);
  info = PLA_Mvector_create_conf_to(b, 1, &c_y); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &h_dot_x); CheckError(info);
  info = PLA_Mvector_create_conf_to(b, 1, &rho); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &nu); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &tau); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &sigma); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &gamma_z); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &gamma_s); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &hat_nu); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &hat_tau); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &delta_x); CheckError(info);
  info = PLA_Mvector_create_conf_to(b, 1, &delta_y); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &delta_s); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &delta_z); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &delta_g); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &delta_t); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &d); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &ones_n); CheckError(info);
  info = PLA_Mvector_create_conf_to(b, 1, &ones_m); CheckError(info);
  info = PLA_Mvector_create_conf_to(c, 1, &max_ratio); CheckError(info);

  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &b_plus_1); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &c_plus_1); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &x_h_x); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &primal_inf); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &dual_inf); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &primal_obj); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &dual_obj); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &mu); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &alfa); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &plus_one); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &minus_one); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &tmp); CheckError(info);
  info = PLA_Mscalar_create(MPI_INT, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &max_idx); CheckError(info);
/*   printf("done.\n"); fflush(stdout); */
  /* Initialize pointers to local buffers. */
  loc.c = NULL;
  info = PLA_Obj_local_buffer(c, (void **) &addr_buf);
  loc.c = (LOQOfloat *) addr_buf;
  loc.b = NULL;
  info = PLA_Obj_local_buffer(b, (void **) &addr_buf);
  loc.b = (LOQOfloat *) addr_buf;
  loc.x = NULL;
  info = PLA_Obj_local_buffer(*x, (void **) &addr_buf);
  loc.x = (LOQOfloat *) addr_buf;
  loc.l = NULL;
  info = PLA_Obj_local_buffer(l, (void **) &addr_buf);
  loc.l = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(u, (void **) &addr_buf);
  loc.u = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(*g, (void **) &addr_buf);
  loc.g = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(*t, (void **) &addr_buf);
  loc.t = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(*s, (void **) &addr_buf);
  loc.s = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(*z, (void **) &addr_buf);
  loc.z = (LOQOfloat *) addr_buf;

  info = PLA_Obj_local_buffer(delta_x, (void **) &addr_buf);
  loc.delta_x = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(delta_g, (void **) &addr_buf);
  loc.delta_g = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(delta_t, (void **) &addr_buf);
  loc.delta_t = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(delta_z, (void **) &addr_buf);
  loc.delta_z = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(delta_s, (void **) &addr_buf);
  loc.delta_s = (LOQOfloat *) addr_buf;

  info = PLA_Obj_local_buffer(sigma, (void **) &addr_buf);
  loc.sigma = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(nu, (void **) &addr_buf);
  loc.nu = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(tau, (void **) &addr_buf);
  loc.tau = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(hat_nu, (void **) &addr_buf);
  loc.hat_nu = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(hat_tau, (void **) &addr_buf);
  loc.hat_tau = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(gamma_z, (void **) &addr_buf);
  loc.gamma_z = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(gamma_s, (void **) &addr_buf);
  loc.gamma_s = (LOQOfloat *) addr_buf;

  info = PLA_Obj_local_buffer(c_plus_1, (void **) &addr_buf);
  loc.c_plus_1 = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(b_plus_1, (void **) &addr_buf);
  loc.b_plus_1 = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(primal_inf, (void **) &addr_buf);
  loc.primal_inf = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(dual_inf, (void **) &addr_buf);
  loc.dual_inf = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(primal_obj, (void **) &addr_buf);
  loc.primal_obj = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(dual_obj, (void **) &addr_buf);
  loc.dual_obj = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(x_h_x, (void **) &addr_buf);
  loc.x_h_x = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(tmp, (void **) &addr_buf);
  loc.tmp = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(mu, (void **) &addr_buf);
  loc.mu = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(alfa, (void **) &addr_buf);
  loc.alfa = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(max_ratio, (void **) &addr_buf);
  loc.max_ratio = (LOQOfloat *) addr_buf;

  info = PLA_Obj_local_buffer(d, (void **) &addr_buf);
  loc.d = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(c_x, (void **) &addr_buf);
  loc.c_x = (LOQOfloat *) addr_buf;
  info = PLA_Obj_local_buffer(c_y, (void **) &addr_buf);
  loc.c_y = (LOQOfloat *) addr_buf;

  /* Determine local size of vectors and stride. */
  info = PLA_Obj_local_length(*x, &local_n); CheckError(info);
  info = PLA_Obj_local_stride(*x, &local_stride); CheckError(info);
  info = PLA_Obj_local_length(b, &local_m); CheckError(info);
  info = PLA_Obj_local_stride(b, &local_stride2); CheckError(info);

  /* Initial settings */
/*   printf("Initial settings ..."); fflush(stdout); */
  info = PLA_Obj_set_to_one(plus_one); CheckError(info);
  info = PLA_Obj_set_to_minus_one(minus_one); CheckError(info);
  info = PLA_Obj_set_to_one(ones_n); CheckError(info);
  info = PLA_Obj_set_to_one(ones_m); CheckError(info);
  info = PLA_Dot(c, c, c_plus_1); CheckError(info);
  info = PLA_Dot(b, b, b_plus_1); CheckError(info);
/*   info = PLA_Dot(c, ones_n, c_plus_1); CheckError(info); */
/*   info = PLA_Dot(b, ones_m, b_plus_1); CheckError(info); */
  *loc.c_plus_1 = sqrt(*loc.c_plus_1);
  *loc.b_plus_1 = sqrt(*loc.b_plus_1);
  *loc.c_plus_1 += 1.0;
  *loc.b_plus_1 += 1.0;
/*   printf("done.\n"); fflush(stdout); */

  /* Get diagonal terms */
/*   printf("Get diagonal terms..."); fflush(stdout); */
  PLA_Get_diagonal(h_x, &diag_h_x);
  PLA_Copy(h_x, h_x_copy);
/*   printf("done.\n"); fflush(stdout); */

/*   printf("Setting initial values..."); */
  /* Starting point */
  if(restart == 1)
    {
      for(i=0; i<local_n; ++i)
	{
	  idx = i*local_stride;
	  loc.g[idx] = max(ABS(loc.x[idx] - loc.l[idx]), bound);
	  loc.t[idx] = max(ABS(loc.u[idx] - loc.x[idx]), bound);
	}
      /* Compute 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);      
      /* Compute sigma <- c + h_dot_x - A^T*y */
      PLA_Copy_vector(h_dot_x, sigma);
      info = PLA_Gemv(PLA_TRANS, minus_one, a, *y, plus_one, sigma);
      CheckError(info);
      info = PLA_Axpy(plus_one, c, sigma); CheckError(info);
      /* MPI Barrier to wait for changes in sigma. */
      info = MPI_Barrier(comm);
      /* Set s and z */
      for(i=0; i<local_n; ++i)
	{
	  idx = i*local_stride;
	  if(loc.sigma[idx] > 0)
	    {
	      loc.s[idx] = bound;
	      loc.z[idx] = loc.sigma[idx] + bound;
	    }
	  else
	    {
	      loc.s[idx] = bound - loc.sigma[idx];
	      loc.z[idx] = bound;
	    }
	}
    } /* if(restart == 1) */
  else
    {
      /* Use default start settings */

      /* Set h_y <- m x m identity matrix */
      info = PLA_Obj_set_to_zero(h_y); CheckError(info);
      PLA_Set_unit_diagonal(&h_y);
      /* Set h_x <- h_x + n x n identity matrix */
      info = PLA_Obj_set_to_one(d); CheckError(info);
      info = PLA_Axpy(plus_one, diag_h_x, d); CheckError(info);
      PLA_Set_diagonal(&h_x, d);
      /* Set c_x <- c and c_y <- b */
      PLA_Copy_vector(c, c_x);
      PLA_Copy_vector(b, c_y);
      /* Solve the reduced KKT system. */
      parallel_solve_reduced(h_x, h_y, a, c_x, c_y, x, y, &y1t, PREDICTOR);