// D. Brugger, december 2006// $Id: loqo_solver.cpp 7 2006-12-16 16:45:32Z beeblbrox $//// Copyright (C) 2006 Dominik Brugger//// This program is free software; you can redistribute it and/or modify// it under the terms of the GNU General Public License as published by// the Free Software Foundation; either version 2 of the License, or// (at your option) any later version.// // This program is distributed in the hope that it will be useful,// but WITHOUT ANY WARRANTY; without even the implied warranty of// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the// GNU General Public License for more details.//// You should have received a copy of the GNU General Public License along// with this program; if not, write to the Free Software Foundation, Inc.,// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.#include "loqo_solver.h"#ifndef CheckError#define CheckError(n) if(n){printf("line %d, file %s\n",__LINE__,__FILE__);}#endifSolver_LOQO::Solver_LOQO(int n, int q, int m){   // Ensure that n,q are even numbers.  this->n = n % 2 == 0 ? n : n-1;   this->n_old = this->n;   this->q = q % 2 == 0 ? q : q-1;   this->m = m;  // Ensure sane q  this->q = this->q > this->n ? this->n : this->q;  this->init_margin = 0.1;  this->init_iter = 500;  this->precision_violations = 0;  this->opt_precision = 1e-8;  NEXT_RAND = 1;}unsigned int Solver_LOQO::next_rand_pos(){  NEXT_RAND = NEXT_RAND*1103515245L + 12345L;  return NEXT_RAND & 0x7fffffff;}void Solver_LOQO::setup_up(int *work_set) {   for(int i=0; i<n; ++i)     up[i] = get_C(work_set[i]);}void Solver_LOQO::allocate_a(){   a = new double[n]; }void Solver_LOQO::allocate_d(){   d = new double[1]; }void Solver_LOQO::allocate_work_space(){   work_space = new double[5*n+1]; }void Solver_LOQO::setup_low(){   for(int i=0; i<n; ++i)     low[i] = 0; }Solver_Parallel_LOQO::Solver_Parallel_LOQO(int n, int q, int m, MPI_Comm comm, int nprows, 					   int npcols, int nb_distr) : Solver_LOQO(n,q,m){  this->comm = comm;  ierr = MPI_Comm_rank(comm, &this->rank); CheckError(ierr);  ierr = MPI_Comm_size(comm, &this->size); CheckError(ierr);  ierr = PLA_Comm_1D_to_2D(comm, nprows, npcols, &comm); CheckError(ierr);  // Initialize PLAPACK  ierr = PLA_Init(comm); CheckError(ierr);  ierr = PLA_Temp_create(nb_distr, 0, &templ); CheckError(ierr);  // Initialize PLAPACK objs  Q_bb_global = NULL; c_global = NULL; up_global = NULL;   low_global = NULL; a_global = NULL; d_global = NULL;  x = NULL; dist = NULL; g = NULL; t = NULL; yy = NULL; z = NULL;  s = NULL; plus_one = NULL; Q_bb_global_view = NULL;  c_global_view = NULL; up_global_view = NULL;  low_global_view = NULL; a_global_view = NULL;  x_view = NULL; dist_view = NULL; g_view = NULL; t_view = NULL;  z_view = NULL; s_view = NULL; x_loc_view = NULL;  x_loc = NULL; yy_loc = NULL; dist_loc = NULL;  dist_loc_view = NULL;}Solver_Parallel_LOQO::~Solver_Parallel_LOQO(){  // Finalize PLAPACK  ierr = PLA_Finalize(); CheckError(ierr);}class Solver_Parallel_LOQO_NU : public Solver_Parallel_LOQO{public:  Solver_Parallel_LOQO_NU(int n, int q, int m, MPI_Comm comm, int nprows, 			  int npcols, int nb_distr, double nu)     : Solver_Parallel_LOQO(n, q, m, comm, nprows, npcols, nb_distr)  {    this->nu = nu;  };private:  double nu;  void setup_a(int *work_set);  void setup_d(int *not_work_set);  int select_working_set(int *work_set, int *not_work_set);//   double calculate_rho() { si->r = dual[1]; return dual[0]; }  double calculate_rho();};double Solver_Parallel_LOQO_NU::calculate_rho(){  printf("Solver_Parallel_LOQO_NU::calculate_rho called!\n");  int nr_free1 = 0,nr_free2 = 0;  double ub1 = INF, ub2 = INF;  double lb1 = -INF, lb2 = -INF;  double sum_free1 = 0, sum_free2 = 0;//   printf("alpha = ");//   for(int i=0; i<l; ++i)//     printf(" %g",alpha[i]);//   printf("\n");  for(int i=0;i<l;i++)    {      if(y[i]==+1)	{	  if(is_lower_bound(i))	    ub1 = min(ub1,G[i]);	  else if(is_upper_bound(i))	    lb1 = max(lb1,G[i]);	  else	    {	      ++nr_free1;	      sum_free1 += G[i];	    }	}      else	{	  if(is_lower_bound(i))	    ub2 = min(ub2,G[i]);	  else if(is_upper_bound(i))	    lb2 = max(lb2,G[i]);	  else	    {	      ++nr_free2;	      sum_free2 += G[i];	    }	}    }  printf("nr_free1 = %d\n", nr_free1);  printf("sum_free1 = %g\n",sum_free1);  printf("nr_free2 = %d\n", nr_free2);  printf("sum_freee = %g\n",sum_free2);  double r1,r2;  if(nr_free1 > 0)    r1 = sum_free1/nr_free1;  else    r1 = (ub1+lb1)/2;	  if(nr_free2 > 0)    r2 = sum_free2/nr_free2;  else    r2 = (ub2+lb2)/2;	  si->r = (r1+r2)/2;  printf("(r1+r2)/2 = %g\n", (r1+r2)/2);  printf("(r1+r2)/2 = %g\n", (r1-r2)/2);  return (r1-r2)/2;}void Solver_Parallel_LOQO_NU::setup_a(int *work_set){  info("Solver_Parallel_LOQO_NU::setup_a called\n");  // Note that a has to be in column major layout.  for(int i=n_low_loc; i<n_up_loc; ++i)    {      a[(i-n_low_loc)*m] = y[work_set[i]];      a[(i-n_low_loc)*m+1] = 1;    }  ierr = PLA_Obj_set_to_zero(a_global_view); CheckError(ierr);  ierr = PLA_API_begin(); CheckError(ierr);  ierr = PLA_Obj_API_open(a_global_view); CheckError(ierr);  ierr = PLA_API_axpy_matrix_to_global(m, local_n, plus_one, (void *)a, m, 				       a_global_view, 0, n_low_loc);  CheckError(ierr);  ierr = PLA_Obj_API_close(a_global_view); CheckError(ierr);  ierr = PLA_API_end(); CheckError(ierr);  info("Solver_Parallel_LOQO_NU::setup_a done\n");}void Solver_Parallel_LOQO_NU::setup_d(int *not_work_set){  info("Solver_Parallel_LOQO_NU::setup_d called\n");  d[0] = 0;  for(int i=0; i<lmn; ++i)    {      if(fabs(alpha[not_work_set[i]]) > TOL_ZERO)	d[0] -= y[not_work_set[i]]*alpha[not_work_set[i]];    }  d[1] = nu*l;  info("Setting d[1] = %g\n", d[1]); info_flush();  for(int i=0; i<lmn; ++i)    d[1] -= alpha[not_work_set[i]];  ierr = PLA_Obj_set_to_zero(d_global); CheckError(ierr);  ierr = PLA_API_begin(); CheckError(ierr);  ierr = PLA_Obj_API_open(d_global); CheckError(ierr);  if(rank == 0)    ierr = PLA_API_axpy_vector_to_global(m, plus_one, (void *) d, 1, d_global, 					 0); CheckError(ierr);  ierr = PLA_Obj_API_close(d_global); CheckError(ierr);  ierr = PLA_API_end(); CheckError(ierr);  info("Solver_Parallel_LOQO_NU::setup_d done\n"); info_flush();}void Solver_Parallel_LOQO::setup_range(int *range_low, int *range_up, 				       int total_sz){  int local_sz = total_sz/size;  int idx_up = local_sz;  int idx_low = 0;  if(total_sz != 0)    {      for(int i=0; i<size-1; ++i)	{	  range_low[i] = idx_low;	  range_up[i] = idx_up;	  idx_low = idx_up;	  idx_up = idx_low + local_sz + 1;	}      range_low[size-1] = idx_low; range_up[size-1]=total_sz;    }  else    {      for(int i=0; i<size; ++i)	{	  range_low[i] = 0;	  range_up[i] = 0;	}    }}int Solver_Parallel_LOQO_NU::select_working_set(int *work_set, 						int *not_work_set){  printf("selecting working set...");  // reset work status  n = n_old;  for(int i=0; i<l; ++i)    {      work_status[i] = WORK_N;    }  double Gmin1 = INF; double Gmin2 = INF;  double Gmax1 = -INF; double Gmax2 = -INF;  int min1 = -1; int min2 = -1;  int max1 = -1; int max2 = -1;  for(int t=0; t<l; ++t)    {      if(y[t] == +1)	{	  if(!is_upper_bound(t))	    {	      if(G[t] < Gmin1)		{		  Gmin1 = G[t];		  min1 = t;		}	    }	  if(!is_lower_bound(t))	    {	      if(G[t] > Gmax1)		{		  Gmax1 = G[t];		  max1 = t;		}	    }	}      else	{	  if(!is_upper_bound(t))	    {	      if(G[t] < Gmin2)		{		  Gmin2 = G[t];		  min2 = t;		}	    }	  if(!is_lower_bound(t))	    {	      if(G[t] > Gmax2)		{		  Gmax2 = G[t];		  max2 = t;		}	    }	}    }  // check for optimality, max. violating pair.  printf("max(Gmax1-Gmin1,Gmax2-Gmin2) = %g < %g\n", 	 max(Gmax1-Gmin1,Gmax2-Gmin2),eps);  if(max(Gmax1-Gmin1,Gmax2-Gmin2) < eps)    return 1;  // Sort gradient  double *Gtmp = new double[l];  int *pidx = new int[l];  for(int i=0; i<l; ++i)    {      Gtmp[i] = G[i];      pidx[i] = i;    }  quick_sort(Gtmp, pidx, 0, l-1);//   printf("pidx = ");//   for(int i=0; i<l; ++i)//     printf(" %d", pidx[i]);//   printf("\n");  int top1=l-1; int top2=l-1;  int bot1=0; int bot2=0;  int count=0;  // Select a full set initially  int nselect = iter == 0 ? n : q;  while(count < nselect)    {      if(top1 > bot1)	{ 	  while(!( (is_free(pidx[top1]) || is_upper_bound(pidx[top1])) 		   && y[pidx[top1]] == +1))	    {	      if(top1 <= bot1) break;	      --top1;	    } 	  while(!( (is_free(pidx[bot1]) || is_lower_bound(pidx[bot1])) 		   && y[pidx[bot1]] == +1))	    {	      if(bot1 >= top1) break;	      ++bot1;	    }	}      if(top2 > bot2)	{	  while(!( (is_free(pidx[top2]) || is_upper_bound(pidx[top2]))		   && y[pidx[top2]] == -1))	    {	      if(top2 <= bot2) break;	      --top2;	    }	  while(!( (is_free(pidx[bot2]) || is_lower_bound(pidx[bot2]))		   && y[pidx[bot2]] == -1))	    {	      if(bot2 >= top2) break;	      ++bot2;	    }	}      if(top1 > bot1 && top2 > bot2)	{	  if(G[pidx[top1]]-G[pidx[bot1]] > G[pidx[top2]]-G[pidx[bot2]])	    {	      work_status[pidx[top1]] = WORK_B;	      work_status[pidx[bot1]] = WORK_B;	      --top1;	      ++bot1;	    }	  else	    {	      work_status[pidx[top2]] = WORK_B;	      work_status[pidx[bot2]] = WORK_B;	      --top2;	      ++bot2;	    }	  count += 2;	}      else if(top1 > bot1)	{	  work_status[pidx[top1]] = WORK_B;	  work_status[pidx[bot1]] = WORK_B;	  --top1;	  ++bot1;	  count += 2;	}      else if(top2 > bot2)	{	  work_status[pidx[top2]] = WORK_B;	  work_status[pidx[bot2]] = WORK_B;	  --top2;	  ++bot2;	  count += 2;	}      else	break;    } // while(count < nselect)  if(count < n)    {      // Compute subset of indices in previous working set       // which were not yet selected      int j=0;      int *work_count_subset = new int[l-count];      int *subset = new int[l-count];      int *psubset = new int[l-count];      for(int i=0; i<l; ++i)	{	  if(work_status[i] == WORK_N && work_count[i] > -1)	    {	      work_count_subset[j] = work_count[i];	      subset[j] = i;	      psubset[j] = j;	      ++j;	    }	}      quick_sort(work_count_subset, psubset, 0, j-1);      // Fill up with j \in B, 0 < alpha[j] < C      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_free(subset[psubset[i]])) 	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // Fill up with j \in B, alpha[j] = 0      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_lower_bound(subset[psubset[i]]))	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // Fill up with j \in B, alpha[j] = C      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_upper_bound(subset[psubset[i]]))	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // clean up      delete[] work_count_subset;      delete[] subset;      delete[] psubset;    } // if(count < n)  // Setup work_set and not_work_set  // update work_count  int nnew=0; int i=0; int j=0; n=0;  for(int t=0; t<l; ++t)    {      if(work_status[t] == WORK_B)	{	  if(work_count[t] == -1)	    ++nnew;	  work_set[i] = t; ++i; ++n;	  ++work_count[t];	}      else	{	  not_work_set[j] = t; ++j;	  work_count[t] = -1;	}    }  // Update q  printf("nnew = %d\n", nnew);  int kin = nnew;  nnew = nnew % 2 == 0 ? nnew : nnew-1;  int L = n/10 % 2 == 0 ? n/10 : (n/10)-1;  q = min(q, max( max( 10, L ), nnew ) );  printf("q = %d\n", q);  printf("n = %d\n",n);  if(kin == 0)    {      // Increase precision of solver.      if(opt_precision > 1e-20)	opt_precision /= 100;      else	{	  info("Error: Unable to select a suitable working set!!!\n");	  return 1; 	}    }  // clean up  delete[] Gtmp;  delete[] pidx;  printf("done.\n");  return 0;  }void Solver_Parallel_LOQO::allocate_a()