// 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);