parallel_loqo.c

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

// D. Brugger, september 2006
// parallel_loqo.c - parallel implementation of LOQO qp solver using PLAPACK.
// $Id: parallel_loqo.c 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 <PLA.h>
#include "PLA_missing_decls.h"
#include "parallel_loqo.h"

#define PREDICTOR 1
#define CORRECTOR 2

#define	max(A, B)	((A) > (B) ? (A) : (B))
#define	ABS(A)  	((A) > 0 ? (A) : (-(A)))
#define sqr(A)          ((A) * (A))

#ifndef CheckError
#define CheckError(n) if(n){printf("line %d, file %s, error %d\n",__LINE__,__FILE__,n); }
#endif

/* Auxiliary functions for setting and getting diagonal values of a 
   matrix.
*/

void PLA_Set_unit_diagonal(PLA_Obj *a)
{
  PLA_Obj a_cur = NULL, a11 = NULL; 
  void *a11_buf;
  int size, info;
  int local_n, local_m;
  info = PLA_Obj_view_all(*a, &a_cur); CheckError(info);
  while(1)
    {
      info = PLA_Obj_global_length(a_cur, &size); CheckError(info);
      if(size == 0) break;
      info = PLA_Obj_split_4(a_cur, 1, 1, &a11, PLA_DUMMY, PLA_DUMMY, &a_cur);
      CheckError(info);
      info = PLA_Obj_local_length(a11, &local_m); CheckError(info);
      info = PLA_Obj_local_width(a11, &local_n); CheckError(info);
      if(local_n == 1 && local_m == 1)
	{
	  info = PLA_Obj_local_buffer(a11, (void **) &a11_buf); 
	  CheckError(info);
	  *((LOQOfloat*) a11_buf) = 1;
	}
    }
}

void print_matrix(PLA_Obj A)
{
  void *local_buf = NULL;
  int local_m, local_n, local_ldim;
  int i,j;
  PLA_Obj_local_length( A, &local_m );
  PLA_Obj_local_width(  A, &local_n );
  PLA_Obj_local_buffer( A, (void **) &local_buf );
  PLA_Obj_local_ldim(   A, &local_ldim );
  for (j=0; j<local_n; j++ )
    {
      for (i=0; i<local_m; i++ )
	printf(" %g", ((LOQOfloat *) local_buf)[ j*local_ldim + i ]);
      printf("\n");
    }
  fflush(stdout);
}

void print_vector(PLA_Obj v)
{
  void *addr_buf = NULL;
  int info;
  int local_stride, local_n, i;
  info = PLA_Obj_local_length(v, &local_n); CheckError(info);
  info = PLA_Obj_local_stride(v, &local_stride); CheckError(info);  
  info = PLA_Obj_local_buffer(v, (void **) &addr_buf); CheckError(info);
  for(i=0; i<local_n; ++i)
    {
      printf(" %g", ((LOQOfloat *) addr_buf)[local_stride*i]);
    }
  printf("\n"); fflush(stdout);
}

/* Assuming vectors have same length + layout */
void PLA_Copy_vector(PLA_Obj from, PLA_Obj to)
{
  PLA_Template templ;
  MPI_Comm comm;
  void *from_buf = NULL, *to_buf = NULL;
  int info;
  int local_stride, local_n, i;
  info = PLA_Obj_template(from, &templ); CheckError(info);
  info = PLA_Temp_comm_all(templ, &comm); CheckError(info);
  info = PLA_Obj_local_length(from, &local_n); CheckError(info);
  info = PLA_Obj_local_stride(from, &local_stride); CheckError(info);  
  info = PLA_Obj_local_buffer(from, (void **) &from_buf); CheckError(info);
  info = PLA_Obj_local_buffer(to, (void **) &to_buf); CheckError(info);
  for(i=0; i<local_n; ++i)
    {
      ((LOQOfloat*) to_buf)[i*local_stride] = ((LOQOfloat *) from_buf)[i*local_stride];
    }
  info = MPI_Barrier(comm); CheckError(info);
}

void PLA_Shift_diagonal(PLA_Obj *a, LOQOfloat shift)
{
  PLA_Obj a_cur = NULL, a11 = NULL; 
  void *a11_buf;
  int size, info;
  int local_m, local_n;
  info = PLA_Obj_view_all(*a, &a_cur);
  while(1)
    {
      info = PLA_Obj_global_length(a_cur, &size);
      if(size == 0) break;
      info = PLA_Obj_split_4(a_cur, 1, 1, &a11, PLA_DUMMY, PLA_DUMMY, &a_cur);
      info = PLA_Obj_local_length(a11, &local_m); CheckError(info);
      info = PLA_Obj_local_width(a11, &local_n); CheckError(info);
      if(local_m == 1 && local_n == 1)
	{
	  info = PLA_Obj_local_buffer(a11, (void **) &a11_buf);
	  *((LOQOfloat*) a11_buf) += shift;
	}
    }
}

void PLA_Set_diagonal(PLA_Obj *a, PLA_Obj d)
{
  PLA_Obj a_cur = NULL, d_cur = NULL, a11 = NULL, d1 = NULL; 
  void *a11_buf; void *d1_buf;
  int size, info;
  int local_m, local_n, local_m2;
  info = PLA_Obj_view_all(*a, &a_cur);
  info = PLA_Obj_view_all(d, &d_cur);
  while(1)
    {
      info = PLA_Obj_global_length(a_cur, &size);
      if(size == 0) break;
      info = PLA_Obj_split_4(a_cur, 1, 1, &a11, PLA_DUMMY, PLA_DUMMY, &a_cur);
      info = PLA_Obj_horz_split_2(d_cur, 1, &d1, &d_cur);
      info = PLA_Obj_local_length(a11, &local_m); CheckError(info);
      info = PLA_Obj_local_width(a11, &local_n); CheckError(info);
      info = PLA_Obj_local_length(d1, &local_m2); CheckError(info);
      if(local_m == 1 && local_n == 1 && local_m2 == 1)
	{
	  info = PLA_Obj_local_buffer(a11, (void **) &a11_buf);
	  info = PLA_Obj_local_buffer(d1, (void **) &d1_buf);
	  *((LOQOfloat*) a11_buf) = *((LOQOfloat*) d1_buf);
	}
    }
}

void PLA_Get_diagonal(PLA_Obj a, PLA_Obj *d)
{
  PLA_Obj a_cur = NULL, d_cur = NULL, a11 = NULL, d1 = NULL; 
  void *a11_buf; void *d1_buf;
  int size, info;
  int local_m, local_n, local_m2;
  info = PLA_Obj_view_all(a, &a_cur);
  info = PLA_Obj_view_all(*d, &d_cur);
  while(1)
    {
      info = PLA_Obj_global_length(a_cur, &size);
      if(size == 0) break;
      info = PLA_Obj_split_4(a_cur, 1, 1, &a11, PLA_DUMMY, PLA_DUMMY, &a_cur);
      info = PLA_Obj_horz_split_2(d_cur, 1, &d1, &d_cur);
      info = PLA_Obj_local_length(a11, &local_m); CheckError(info);
      info = PLA_Obj_local_width(a11, &local_n); CheckError(info);
      info = PLA_Obj_local_length(d1, &local_m2); CheckError(info);
      if(local_m == 1 && local_n == 1 && local_m2 == 1)
	{
	  info = PLA_Obj_local_buffer(a11, (void **) &a11_buf);
	  info = PLA_Obj_local_buffer(d1, (void **) &d1_buf);
	  *((LOQOfloat*) d1_buf) = *((LOQOfloat*) a11_buf);
	}
    }
}

/* Cholesky decomposition with manteuffel shifting. */
int PLA_Choldc(PLA_Obj a)
{
  int shifted, i, idx;
  int info = 0;
  LOQOfloat rs;
  LOQOfloat shift_amount = 0.0;
  LOQOfloat one = 1.0;
  LOQOfloat *array = NULL;
  int error,parameters,sequential,r12;
  do{
    shifted = 0;
    /* This is necessary as PLA_Chol, returns the index
       of a non-positive pivot in info only if error checking
       is disabled. */
    info = PLA_Get_error_checking(&error, &parameters, &sequential, &r12);
    info = PLA_Set_error_checking(FALSE, FALSE, FALSE, FALSE);
    info = PLA_Chol(PLA_LOWER_TRIANGULAR, a);
    if(info != 0){
      /* We need a zero based index */
      idx = info-1;
      info = PLA_API_begin(); CheckError(info);
      info = PLA_Obj_API_open(a); CheckError(info);
      array = (LOQOfloat*) malloc((idx+1)*sizeof(LOQOfloat));
      /* Set array to zero, because of axpy op. */
      for(i=0; i<idx+1; ++i)
	array[i] = 0;
      info = PLA_API_axpy_global_to_matrix(1,idx+1,&one,a,idx,0,
					   (void *)array,1); CheckError(info);
      info = PLA_Obj_API_close(a); CheckError(info);
      info = PLA_API_end(); CheckError(info);
      /* Compute part of row sum. */
      rs = 0.0;
      for(i=0; i<idx; ++i) /* Add pivot value below */
	rs += fabs(array[i]);
      free(array);

      info = PLA_API_begin(); CheckError(info);
      info = PLA_Obj_API_open(a); CheckError(info);
      array = (LOQOfloat*) malloc((idx+1)*sizeof(LOQOfloat));
      /* Set array to zero, because of axpy op. */
      for(i=0; i<idx+1; ++i)
	array[i] = 0;
      info = PLA_API_axpy_global_to_matrix(idx+1,1,&one,a,idx,idx,
					   (void *)array,1); CheckError(info);
      info = PLA_Obj_API_close(a); CheckError(info);
      info = PLA_API_end(); CheckError(info);
      /* Complete computation of row sum. */
      for(i=0; i<idx+1; ++i)
	rs += fabs(array[i]);
      shift_amount = max(rs,1.1*shift_amount);
      printf("using shift_amount = %g\n", shift_amount); fflush(stdout);
      free(array);

      PLA_Shift_diagonal(&a, shift_amount);
      shifted = 1;
    }
    /* Restore error checking state. */
    info = PLA_Set_error_checking(error, parameters, sequential, r12);
  } while(shifted);

  return info;
}

/* Solves reduced KKT system:
   | -Q1 A^T | | delta_alpha | = | c1 |
   | A   Q2  | | delta_h     |   | c2 |
   Dimension of Q1 is n x n. Dimension of A is m x n.
   Y1T and stores intermediary result computed during predictor step.
   Dimension of Y1T is m x n.
*/
void parallel_solve_reduced(PLA_Obj Q1, PLA_Obj Q2, PLA_Obj A, PLA_Obj c1, PLA_Obj c2, 
			    PLA_Obj *delta_alpha, PLA_Obj *delta_h, PLA_Obj *Y1T,
			    int step)
{
  int info;
  PLA_Obj Y2 = NULL, minus_one = NULL, plus_one = NULL;
  PLA_Template templ = NULL;
  info = PLA_Obj_template(Q1, &templ); CheckError(info);
  info = PLA_Mvector_create_conf_to(*delta_alpha, 1, &Y2); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &plus_one); CheckError(info);
  info = PLA_Obj_set_to_one(plus_one); CheckError(info);
  info = PLA_Mscalar_create(MPIfloat, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, 
			    templ, &minus_one); CheckError(info);
  info = PLA_Obj_set_to_minus_one(minus_one); CheckError(info);
  if(step == PREDICTOR)
    {
      /* Compute cholesky decomposition: Q1 <- L1*L1^T
         overwriting lower triangular part of Q1 */
/*       info = PLA_Chol(PLA_LOWER_TRIANGULAR, Q1); CheckError(info); */
      info = PLA_Choldc(Q1); CheckError(info);
      /* Compute: Y1^T <- A * (L1^-1)^T.
	 Dimension of Y1^T is m x n. */
      info = PLA_Copy(A, *Y1T); CheckError(info);
      info = PLA_Trsm(PLA_SIDE_RIGHT, PLA_LOWER_TRIANGULAR, PLA_TRANS,
		      PLA_NONUNIT_DIAG, plus_one, Q1, *Y1T); CheckError(info);
      /* Compute cholesky decomposition: (Q2 + Y1^T*Y1) <- L2*L2^T. */
      info = PLA_Syrk(PLA_LOWER_TRIANGULAR, PLA_NO_TRANS, plus_one, *Y1T,
		      plus_one, Q2); CheckError(info);
/*       info = PLA_Chol(PLA_LOWER_TRIANGULAR, Q2); CheckError(info);  */
      info = PLA_Choldc(Q2); CheckError(info);
    }
  /* Compute Y2 <- L1^-1*c1. 
     Dimension of Y2 is n x 1. */
  PLA_Copy_vector(c1, Y2);
  info = PLA_Trsv(PLA_LOWER_TRIANGULAR, PLA_NO_TRANS, PLA_NONUNIT_DIAG,
		  Q1, Y2); CheckError(info);
  /* Compute delta_h <- c2 + Y1^T*Y2. */
  PLA_Copy_vector(c2, *delta_h);
  info = PLA_Gemv(PLA_NO_TRANS, plus_one, *Y1T, Y2, plus_one, *delta_h);
  CheckError(info);
  /* Compute delta_h <- L2^-1*(c2 + Y1^T*Y2). */
  info = PLA_Trsv(PLA_LOWER_TRIANGULAR, PLA_NO_TRANS, PLA_NONUNIT_DIAG,
		  Q2, *delta_h); CheckError(info);
  /* Finally compute delta_h <- L2^-T * (L2^-1*(c2 + Y1^T*Y2). */
  info = PLA_Trsv(PLA_LOWER_TRIANGULAR, PLA_TRANS, PLA_NONUNIT_DIAG,
		  Q2, *delta_h); CheckError(info);
  /* Compute delta_alpha <- Y1*delta_h - Y2. */
  PLA_Copy_vector(Y2, *delta_alpha);
  info = PLA_Gemv(PLA_TRANS, plus_one, *Y1T, *delta_h, minus_one, 
		  *delta_alpha); CheckError(info);
  /* Compute delta_alpha <- L1^-T * (Y1*delta_h - Y2). */
  info = PLA_Trsv(PLA_LOWER_TRIANGULAR, PLA_TRANS, PLA_NONUNIT_DIAG, Q1,
		  *delta_alpha); CheckError(info);
  info = PLA_Obj_free(&Y2); CheckError(info);
  info = PLA_Obj_free(&plus_one); CheckError(info);
  info = PLA_Obj_free(&minus_one); CheckError(info);