void DistMat ( double *x, double *y, int n,
	       double *dist ) {
  int i, j;

  for ( i = 0; i < n; i++ ) {
    for ( j = 0; j < n; j++ ) {
      dist[j * n + i] = sqrt ( pow( x[i] - x[j], 2.0) + 
			       pow( y[i] - y[j], 2.0) );
    }
  }
}


double Timing()
{
    struct timeval t;
    gettimeofday(&t, NULL);
    return ((double)t.tv_sec + (double)t.tv_usec * 1e-6);
}

void PrintDoubleArrayAsMatrix ( double * array, int m, int n )
{
    int i, j;

    for ( i=0; i<m; i++)
    {
    	for ( j=0; j<n; j++)
    		printf ( "%.3lf ", array[j*n+i] );
        printf ( "\n" );
    }
}

void WriteArrayToFile ( double * array, int m, int n, char* fileName) {
  int i, j;
  FILE *writeFile;

  double *ptr = array;
    
  writeFile = fopen ( fileName, "w" );
      
  fprintf ( writeFile, "%d %d\n", m, n );

  for ( i=0; i<m; i++) {
    for ( j=0; j<n; j++) {
      fprintf ( writeFile, "%lf ", *ptr );
      fprintf ( writeFile, "\n" );
      ptr++;;
    }
  }

  fclose ( writeFile );
  //free ( writeFile );
}

void DataPrep (int me, char *fileName,
	       PLA_Obj *X, PLA_Obj *Z, PLA_Obj *Dist) {
  /*
    The data in fileName is stored this way: z, xdim, ydim
  */
  int nrow, ncol, n, p;
  int i, j, zero_or_one, local_stride, local_ldim;
  double *array, *distArr, dOne = 1.0;
  FILE *file;
  PLA_Obj X1 = NULL, X2 = NULL;

  if ( 0 == me ) {
    file = fopen(fileName, "r");
    fscanf ( file, "%d %d", &nrow, &ncol );
    n = nrow;
    p = ncol;
    array   = (double *) calloc (n * p, sizeof(double));
    distArr = (double *) calloc (n * n, sizeof(double));
    for ( i = 0; i < nrow; i++)
      for ( j = 0; j < ncol; j++)
	fscanf ( file, "%lf", &array[j * nrow + i] ); //important! 
    
    DistMat ( array + n, array + 2 * n, n, distArr );

    fclose( file );
  }
  MPI_Bcast( &p, 1, MPI_INT, 0, MPI_COMM_WORLD );
  MPI_Bcast( &n, 1, MPI_INT, 0, MPI_COMM_WORLD );

  PLA_Temp_zero_or_one( templat, &zero_or_one );

  /* create plapack matrix objects */
  PLA_Matrix_create  ( MPI_DOUBLE, n, n, templat, 0, 0, Dist );
  PLA_Mvector_create ( MPI_DOUBLE, n, 1, templat, 0, Z );
  PLA_Matrix_create  ( MPI_DOUBLE, n, p, templat, 0, 0, X );
  
  PLA_Obj_global_width ( *X, &ncol);

  PLA_Obj_vert_split_2 ( *X, 1, &X1, &X2 );
  
  PLA_Obj_set_to_zero ( *Z );
  PLA_Obj_set_to_one  ( X1 ); //intercept in the design matrix
  PLA_Obj_set_to_zero ( X2 );
  PLA_Obj_set_to_zero ( *Dist );

  // load value to matrix
  PLA_API_begin();
  PLA_Obj_API_open ( *Z ); PLA_Obj_API_open (X2); PLA_Obj_API_open (*Dist);

  if ( 0 == me ) {
    local_stride = 1; //correct!
    PLA_API_axpy_vector_to_global ( n, &dOne, array, local_stride, *Z, zero_or_one);
    local_ldim = n;
    PLA_API_axpy_matrix_to_global ( n, p - 1, &dOne, array + local_ldim, local_ldim, X2, zero_or_one, zero_or_one );
    PLA_API_axpy_matrix_to_global ( n, n, &dOne, distArr, local_ldim, *Dist, zero_or_one, zero_or_one );

    free ( array ); free ( distArr );
  }
    
  PLA_Obj_API_close ( *Z ); PLA_Obj_API_close (X2); PLA_Obj_API_close (*Dist);
  PLA_API_end();

  PLA_Obj_free ( &X1 ); PLA_Obj_free ( &X2 );
}

void ReadFromFile(int me, PLA_Obj *obj, char *fileName, int objtype) {
    int         m, n, i, j, zero_or_one;
    int         local_stride, local_ldim;
    double       *array, f_one = 1.0;
    FILE        *file;
    double	start, end;
    // PLA_Template                    templat = NULL;

    //PLA_Obj_template ( *obj, &templat );
    /* root processor only, load data from file into a temp array */
    if (0 == me) 
    {
      printf ( "in readfromfile\n");
        file = fopen(fileName, "r");
	if (objtype == PLA_MVECTOR) {
	  fscanf( file, "%d", &m);
	  n = 1;
	}
	else {
	  fscanf ( file, "%d %d", &m, &n );
	}
    	array = (double *) calloc (m * n, sizeof(double));

        /* load file into temp array */
    	for ( i=0; i<m; i++)
    	for ( j=0; j<n; j++)
    	    fscanf ( file, "%lf", &array[i*n+j] );

    	fclose( file );
    }

    MPI_Bcast( &m, 1, MPI_INT, 0, MPI_COMM_WORLD );
    MPI_Bcast( &n, 1, MPI_INT, 0, MPI_COMM_WORLD );

    /* get indexing (zero_or_one) from template */
    PLA_Temp_zero_or_one( templat, &zero_or_one );

    /* create plapack matrix objects */
    if (objtype == PLA_MVECTOR) {
      PLA_Mvector_create(MPI_DOUBLE, m, n, templat, 0, obj);
    }
    else {
      PLA_Matrix_create( MPI_DOUBLE, m, n, templat, 0, 0, obj);
    }

    PLA_Obj_set_to_zero(*obj);

    // load value to matrix
    PLA_API_begin();
    PLA_Obj_API_open ( *obj );

    if ( 0 == me )
      {
	if (objtype == PLA_MVECTOR) {
	  local_stride = 1; //correct!
	  PLA_API_axpy_vector_to_global ( m, &f_one, array, local_stride, *obj, zero_or_one);
	}
	else {
	  local_ldim = m;
	  PLA_API_axpy_matrix_to_global ( m, n, &f_one, array, local_ldim, *obj, zero_or_one, zero_or_one );
	}
	free ( array );
      }
    
    PLA_Obj_API_close ( *obj );
    PLA_API_end();
}

typedef double fun(double, double);

double exponential(double dist, double phi) {
  return exp(- phi * dist) ;
}

double myLog(double a, double phi) {
  return log(a);
}

void myPLA_local_apply ( PLA_Obj matrix, fun *f, double phi)
{
    int local_length, local_width, local_ldim, local_stride, i, j;
    int objtype, jump = 0;
    void *local_buffer;
    double *buf;

    MPI_Datatype datatype;

    //if (me == 0)    printf("in myPLA_local_apply\n");
    /* extract  matrix information */
    PLA_Obj_datatype( matrix, &datatype);
    PLA_Obj_objtype( matrix, &objtype);
    PLA_Obj_local_length( matrix, &local_length );
    PLA_Obj_local_width( matrix, &local_width );
    PLA_Obj_local_buffer( matrix, &local_buffer );
    PLA_Obj_local_ldim( matrix, &local_ldim );
    PLA_Obj_local_stride( matrix, &local_stride);

    buf = (double *) local_buffer;
    if (objtype == PLA_MATRIX) jump = local_ldim;
    else jump == local_stride;
    /* apply function elment-wise */
    
    for (i = 0; i < local_length; i++) {
      for (j = 0; j < local_width; j++) {
	//Please, storage is colum major!
	buf[j * jump + i] = f( buf[j * jump + i], phi);
      }
    }
}

double sum_vector(PLA_Obj v) {
  PLA_Obj dOnes = NULL, mean = NULL;
  int n, dummy;
  double mean_val;
  //PLA_Template                    templat = NULL;

  //PLA_Obj_template ( v, &templat );

  PLA_Obj_global_length(v, &n);
  PLA_Mvector_create_conf_to(v, 1, &dOnes);
  PLA_Obj_set_to_one(dOnes);


  PLA_Mscalar_create(MPI_DOUBLE, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, templat, &mean);
  PLA_Dot(v, dOnes, mean);

  PLA_Obj_get_local_contents(mean, PLA_NO_TRANS, &dummy, &dummy, &mean_val, 1, 1);
  //mean_val = mean_val / (double) n;

  PLA_Obj_free(&dOnes);  PLA_Obj_free(&mean); 

  return mean_val;
}

double mean_vector(PLA_Obj v) {
  int n;

  PLA_Obj_global_length(v, &n);
  return sum_vector(v) / (double) n;
}

double var_vector(PLA_Obj v) {
  PLA_Obj minus_one = NULL, one = NULL;
  PLA_Obj alpha = NULL, v_cpy = NULL, var = NULL;
  int n, dummy;
  double minus_mean = - mean_vector(v);
  double var_val;
  //PLA_Template                    templat = NULL;

  //PLA_Obj_template ( v, &templat );

  PLA_Obj_global_length(v, &n);
  PLA_Mvector_create_conf_to(v, 1, &v_cpy);
  //This is not implemented in PLAPACK!!!
  //PLA_Copy(v, v_cpy);
  PLA_Obj_set(v_cpy, MPI_DOUBLE, &minus_mean);
  //PLA_Global_show("v_cpy", v_cpy, "%lf", "end of v_cpy");
  PLA_Create_constants_conf_to(v, &minus_one, NULL, &one);

  PLA_Axpy(one, v, v_cpy);
  PLA_Mscalar_create(MPI_DOUBLE, PLA_ALL_ROWS, PLA_ALL_COLS, 1, 1, templat, &var);
  PLA_Dot(v_cpy, v_cpy, var);
  PLA_Obj_get_local_contents(var, PLA_NO_TRANS, &dummy, &dummy, &var_val, 1, 1);
  var_val = var_val / (double) (n - 1);

  PLA_Obj_free(&minus_one);   PLA_Obj_free(&alpha);
  PLA_Obj_free(&v_cpy);   PLA_Obj_free(&var);

  return var_val;
}

double logProd_diag_matrix(PLA_Obj matrix) {
  /*
    This function computes \sum log(m_ii)
  */
  PLA_Obj diag = NULL;
  double rval;

  PLA_Mvector_create_conf_to(matrix, 1, &diag);
  PLA_Copy_diag_to_mv(matrix, diag);
  myPLA_local_apply(diag, myLog, 0.0);
  rval = sum_vector(diag);
  PLA_Obj_free(&diag);

  return rval;
}

double value_of_scalar(PLA_Obj mscalar) {
  void *buffer;
  double *bufp;
  PLA_Obj_local_buffer(mscalar, &buffer);
  bufp = (double *) buffer;
  return *bufp;
}

void myPLA_copy_vector_to_local(PLA_Obj v, double *loc) {
  int size, local_stride = 1;
  double dOne = 1.0;

  int i;

  PLA_Obj_global_length(v, &size);

  for (i = 0; i < size; i++) loc[i] = 0.0;

  PLA_API_begin();
  PLA_Obj_API_open(v);
  PLA_API_axpy_global_to_vector(size, &dOne, v, 0, loc, local_stride);
  PLA_Obj_API_close(v);
  PLA_API_end();
}

void myPLA_copy_to_local (PLA_Obj v, double *loc ) {
  int size, ldim;
  int nrow, ncol;
  PLA_Obj v_msc = NULL;

  PLA_Obj_global_length ( v, &size );
  PLA_Mscalar_create_conf_to ( v, PLA_ALL_ROWS, PLA_ALL_COLS, &v_msc );
  PLA_Copy ( v, v_msc );
  PLA_Obj_local_ldim ( v_msc, &ldim);

  PLA_Obj_get_local_contents ( v_msc, PLA_NO_TRANS, &nrow, &ncol, loc, ldim, 1 );
}


void myPLA_copy_vector_to_global(PLA_Obj v, double *loc) {
  int size, local_stride = 1;
  double dOne = 1.0;
  int me;

  PLA_Obj_global_length(v, &size);
  MPI_Comm_rank(MPI_COMM_WORLD, &me);

  PLA_Obj_set_to_zero(v);

  PLA_API_begin();
  PLA_Obj_API_open(v);
  if (me == 0) {
    PLA_API_axpy_vector_to_global(size, &dOne, loc, local_stride, v, 0);
  }
  PLA_Obj_API_close(v);
  PLA_API_end();
}

void myPLA_copy_matrix_to_local(PLA_Obj matrix, double *loc) {
  int size_row, size_col, local_ldim;
  double dOne = 1.0;
  int i;
  PLA_Obj_global_length(matrix, &size_row);
  PLA_Obj_global_width(matrix, &size_col);
  local_ldim = size_row;

  for (i = 0; i < size_row * size_col; i++) loc[i] = 0.0;
  PLA_API_begin();
  PLA_Obj_API_open(matrix);
  PLA_API_axpy_global_to_matrix(size_row, size_col, &dOne, matrix, 0, 0, loc, local_ldim);
  PLA_Obj_API_close(matrix);
  PLA_API_end();
}

typedef double fun_ptr(double, double *);

double my_exponential(double dist, double *par) {
  double phi = par[0], S = par[1];
  if (dist == 0) return 1.0;
  else return (1.0 - S) * exp(- phi * dist) ;
}

void myPLA_local_apply_ptr ( PLA_Obj matrix, fun_ptr *f, double *par)
{
    int local_length, local_width, local_ldim, local_stride, i, j;
    int objtype, jump = 0;
    void *local_buffer;
    double *buf;

    MPI_Datatype datatype;

    PLA_Obj_datatype( matrix, &datatype);
    PLA_Obj_objtype( matrix, &objtype);
    PLA_Obj_local_length( matrix, &local_length );
    PLA_Obj_local_width( matrix, &local_width );
    PLA_Obj_local_buffer( matrix, &local_buffer );
    PLA_Obj_local_ldim( matrix, &local_ldim );
    PLA_Obj_local_stride( matrix, &local_stride);

    buf = (double *) local_buffer;
    if (objtype == PLA_MATRIX) jump = local_ldim;
    else jump == local_stride;
    
    for (i = 0; i < local_length; i++) {
      for (j = 0; j < local_width; j++) {
	//Please, storage is colum major!
	buf[j * jump + i] = f( buf[j * jump + i], par);
      }
    }
}