text:
/***************************************************************************/
/*                                                                         */
/* matrix.c - Matrix class for mruby                                       */
/* Copyright (C) 2015 Paolo Bosetti                                        */
/* paolo[dot]bosetti[at]unitn.it                                           */
/* Department of Industrial Engineering, University of Trento              */
/*                                                                         */
/* This library is free software.  You can redistribute it and/or          */
/* modify it under the terms of the GNU GENERAL PUBLIC LICENSE 2.0.        */
/*                                                                         */
/* This library 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           */
/* Artistic License 2.0 for more details.                                  */
/*                                                                         */
/* See the file LICENSE                                                    */
/*                                                                         */
/***************************************************************************/

#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include "matrix.h"
#include "vector.h"

#pragma mark -
#pragma mark • Utilities

// Garbage collector handler, for play_data struct
// if play_data contains other dynamic data, free it too!
// Check it with GC.start
void matrix_destructor(mrb_state *mrb, void *p_) {
  gsl_matrix *v = (gsl_matrix *)p_;
  gsl_matrix_free(v);
};

// Creating data type and reference for GC, in a const struct
const struct mrb_data_type matrix_data_type = {"matrix_data",
                                               matrix_destructor};

// Utility function for getting the struct out of the wrapping IV @data
void mrb_matrix_get_data(mrb_state *mrb, mrb_value self, gsl_matrix **data) {
  mrb_value data_value;
  data_value = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@data"));

  // Loading data from data_value into p_data:
  Data_Get_Struct(mrb, data_value, &matrix_data_type, *data);
  if (!*data)
    mrb_raise(mrb, E_RUNTIME_ERROR, "Could not access @data");
}

#pragma mark -
#pragma mark • Initializations and setup

// Data Initializer C function (not exposed!)
static void mrb_matrix_init(mrb_state *mrb, mrb_value self, mrb_int n,
                            mrb_int m) {
  mrb_value data_value; // this IV holds the data
  gsl_matrix *p_data;   // pointer to the C struct

  data_value = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@data"));

  // if @data already exists, free its content:
  if (!mrb_nil_p(data_value)) {
    Data_Get_Struct(mrb, data_value, &matrix_data_type, p_data);
    free(p_data);
  }
  // Allocate and zero-out the data struct:
  p_data = gsl_matrix_calloc(n, m);
  if (!p_data)
    mrb_raise(mrb, E_RUNTIME_ERROR, "Could not allocate @data");

  // Wrap struct into @data:
  mrb_iv_set(
      mrb, self, mrb_intern_lit(mrb, "@data"), // set @data
      mrb_obj_value(                           // with value hold in struct
          Data_Wrap_Struct(mrb, mrb->object_class, &matrix_data_type, p_data)));
}

static mrb_value mrb_matrix_initialize(mrb_state *mrb, mrb_value self) {
  mrb_int n, m;
  mrb_get_args(mrb, "ii", &n, &m);

  // Call strcut initializer:
  mrb_matrix_init(mrb, self, n, m);
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@nrows"), mrb_fixnum_value(n));
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@ncols"), mrb_fixnum_value(m));
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@format"),
             mrb_str_new_cstr(mrb, "%10.3f"));
  return mrb_nil_value();
}

static mrb_value mrb_matrix_dup(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_matrix *p_mat = NULL, *p_mat_other = NULL;
  mrb_value args[2];

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  args[0] = mrb_fixnum_value(p_mat->size1);
  args[1] = mrb_fixnum_value(p_mat->size2);
  other = mrb_obj_new(mrb, mrb_class_get(mrb, "Matrix"), 2, args);
  mrb_matrix_get_data(mrb, other, &p_mat_other);
  gsl_matrix_memcpy(p_mat_other, p_mat);
  return other;
}

static mrb_value mrb_matrix_rnd_fill(mrb_state *mrb, mrb_value self) {
  gsl_matrix *p_mat = NULL;
  const gsl_rng_type *T;
  gsl_rng *r;
  mrb_int h, k;

  mrb_matrix_get_data(mrb, self, &p_mat);

  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  for (k = 0; k < p_mat->size2; k++) {
    for (h = 0; h < p_mat->size1; h++) {
      gsl_matrix_set(p_mat, h, k, gsl_rng_uniform(r));
    }
  }
  return self;
}

static mrb_value mrb_matrix_all(mrb_state *mrb, mrb_value self) {
  mrb_float v;
  gsl_matrix *p_mat = NULL;

  mrb_get_args(mrb, "f", &v);
  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  gsl_matrix_set_all(p_mat, v);
  return self;
}

static mrb_value mrb_matrix_zero(mrb_state *mrb, mrb_value self) {
  gsl_matrix *p_mat = NULL;

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  gsl_matrix_set_zero(p_mat);
  return self;
}

static mrb_value mrb_matrix_identity(mrb_state *mrb, mrb_value self) {
  gsl_matrix *p_mat = NULL;

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  gsl_matrix_set_identity(p_mat);
  return self;
}

#pragma mark -
#pragma mark • Tests

static mrb_value mrb_matrix_equal(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_matrix *p_mat, *p_mat_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  mrb_matrix_get_data(mrb, other, &p_mat_other);
  if (1 == gsl_matrix_equal(p_mat, p_mat_other))
    return mrb_true_value();
  else
    return mrb_false_value();
}

#pragma mark -
#pragma mark • Accessors

static mrb_value mrb_matrix_get_row(mrb_state *mrb, mrb_value self) {
  mrb_int i, n;
  mrb_value result, args[1];
  gsl_matrix *p_mat = NULL;
  gsl_vector *p_vec = NULL;

  n = mrb_get_args(mrb, "|i", &i);
  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);

  if (n == 1) {
    if (i >= p_mat->size1) {
      mrb_raise(mrb, E_MATRIX_ERROR, "matrix index out of range!");
    }
    args[0] = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@ncols"));
    result = mrb_obj_new(mrb, mrb_class_get(mrb, "Vector"), 1, args);
    mrb_vector_get_data(mrb, result, &p_vec);
    gsl_matrix_get_row(p_vec, p_mat, i);
  }
  else {
    result = mrb_ary_new_capa(mrb, p_mat->size1);
    for (i = 0; i < p_mat->size1; i++) {
      mrb_ary_push(mrb, result, mrb_funcall(mrb, self, "get_row", 1, mrb_fixnum_value(i)));
    }
  }
  return result;
}

static mrb_value mrb_matrix_get_col(mrb_state *mrb, mrb_value self) {
  mrb_int i;
  mrb_value res, args[1];
  gsl_matrix *p_mat = NULL;
  gsl_vector *p_vec = NULL;

  mrb_get_args(mrb, "i", &i);
  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (i >= p_mat->size2) {
    mrb_raise(mrb, E_MATRIX_ERROR, "matrix index out of range!");
  }
  args[0] = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@nrows"));
  res = mrb_obj_new(mrb, mrb_class_get(mrb, "Vector"), 1, args);
  mrb_vector_get_data(mrb, res, &p_vec);
  gsl_matrix_get_col(p_vec, p_mat, i);
  return res;
}

static mrb_value mrb_matrix_get_ij(mrb_state *mrb, mrb_value self) {
  mrb_value result;
  mrb_int i, j;
  gsl_matrix *p_mat = NULL;
  mrb_int n;

  n = mrb_get_args(mrb, "|ii", &i, &j);
  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (n == 2) {
    if (i >= p_mat->size1 || j >= p_mat->size2) {
      mrb_raise(mrb, E_MATRIX_ERROR, "matrix index out of range!");
    }
    result = mrb_float_value(mrb, gsl_matrix_get(p_mat, (size_t)i, (size_t)j));
  } else {
    result = mrb_matrix_get_row(mrb, self);
  }
  return result;
}

static mrb_value mrb_matrix_set_ij(mrb_state *mrb, mrb_value self) {
  mrb_int i, j;
  mrb_float f;
  gsl_matrix *p_mat = NULL;

  mrb_get_args(mrb, "iif", &i, &j, &f);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (i >= p_mat->size1 || j >= p_mat->size2) {
    mrb_raise(mrb, E_MATRIX_ERROR, "matrix index out of range!");
  }
  gsl_matrix_set(p_mat, (size_t)i, (size_t)j, (double)f);
  return mrb_float_value(mrb, f);
}

static mrb_value mrb_matrix_set_row(mrb_state *mrb, mrb_value self) {
  mrb_int i;
  mrb_value other;
  gsl_matrix *p_mat = NULL;
  gsl_vector *p_vec = NULL;

  mrb_get_args(mrb, "io", &i, &other);
  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (i >= p_mat->size1) {
    mrb_raise(mrb, E_MATRIX_ERROR, "matrix row index out of range!");
  }
  mrb_vector_get_data(mrb, other, &p_vec);
  if (p_mat->size2 != p_vec->size) {
    mrb_raise(mrb, E_MATRIX_ERROR, "Size mismatch!");
  }
  gsl_matrix_set_row(p_mat, i, p_vec);
  return self;
}

static mrb_value mrb_matrix_set_col(mrb_state *mrb, mrb_value self) {
  mrb_int i;
  mrb_value other;
  gsl_matrix *p_mat = NULL;
  gsl_vector *p_vec = NULL;

  mrb_get_args(mrb, "io", &i, &other);
  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (i >= p_mat->size2) {
    mrb_raise(mrb, E_MATRIX_ERROR, "matrix col index out of range!");
  }
  mrb_vector_get_data(mrb, other, &p_vec);
  if (p_mat->size1 != p_vec->size) {
    mrb_raise(mrb, E_MATRIX_ERROR, "Size mismatch!");
  }
  gsl_matrix_set_col(p_mat, i, p_vec);
  return self;
}

#pragma mark -
#pragma mark • Properties

static mrb_value mrb_matrix_max(mrb_state *mrb, mrb_value self) {
  gsl_matrix *p_mat = NULL;
  mrb_matrix_get_data(mrb, self, &p_mat);
  return mrb_float_value(mrb, gsl_matrix_max(p_mat));
}

static mrb_value mrb_matrix_min(mrb_state *mrb, mrb_value self) {
  gsl_matrix *p_mat = NULL;
  mrb_matrix_get_data(mrb, self, &p_mat);
  return mrb_float_value(mrb, gsl_matrix_min(p_mat));
}

static mrb_value mrb_matrix_max_index(mrb_state *mrb, mrb_value self) {
  size_t i, j;
  mrb_value res = mrb_ary_new_capa(mrb, 2);
  gsl_matrix *p_mat = NULL;
  mrb_matrix_get_data(mrb, self, &p_mat);
  gsl_matrix_max_index(p_mat, &i, &j);
  mrb_ary_push(mrb, res, mrb_fixnum_value(i));
  mrb_ary_push(mrb, res, mrb_fixnum_value(j));
  return res;
}

static mrb_value mrb_matrix_min_index(mrb_state *mrb, mrb_value self) {
  size_t i, j;
  mrb_value res = mrb_ary_new_capa(mrb, 2);
  gsl_matrix *p_mat = NULL;
  mrb_matrix_get_data(mrb, self, &p_mat);
  gsl_matrix_min_index(p_mat, &i, &j);
  mrb_ary_push(mrb, res, mrb_fixnum_value(i));
  mrb_ary_push(mrb, res, mrb_fixnum_value(j));
  return res;
}

#pragma mark -
#pragma mark • Operations

static mrb_value mrb_matrix_add(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_matrix *p_mat, *p_mat_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Matrix"))) {
    mrb_matrix_get_data(mrb, other, &p_mat_other);
    if (p_mat->size1 != p_mat_other->size1 ||
        p_mat->size2 != p_mat_other->size2) {
      mrb_raise(mrb, E_MATRIX_ERROR, "matrix dimensions don't match!");
    }
    gsl_matrix_add(p_mat, p_mat_other);
  } else if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Numeric"))) {
    gsl_matrix_add_constant(p_mat, mrb_to_flo(mrb, other));
  }
  return self;
}

static mrb_value mrb_matrix_sub(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_matrix *p_mat, *p_mat_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  mrb_matrix_get_data(mrb, other, &p_mat_other);
  if (p_mat->size1 != p_mat_other->size1 ||
      p_mat->size2 != p_mat_other->size2) {
    mrb_raise(mrb, E_MATRIX_ERROR, "matrix dimensions don't match!");
  }
  gsl_matrix_sub(p_mat, p_mat_other);
  return self;
}

static mrb_value mrb_matrix_mul(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_matrix *p_mat, *p_mat_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Matrix"))) {
    mrb_matrix_get_data(mrb, other, &p_mat_other);
    if (p_mat->size1 != p_mat_other->size1 ||
        p_mat->size2 != p_mat_other->size2) {
      mrb_raise(mrb, E_MATRIX_ERROR, "matrix dimensions don't match!");
    }
    gsl_matrix_mul_elements(p_mat, p_mat_other);
  } else if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Numeric"))) {
    gsl_matrix_scale(p_mat, mrb_to_flo(mrb, other));
  }
  return self;
}

static mrb_value mrb_matrix_prod(mrb_state *mrb, mrb_value self) {
  mrb_value other, res;
  gsl_matrix *p_mat, *p_mat_other, *p_mat_res;
  gsl_vector *p_vec_other, *p_vec_res;
  mrb_value args[2];
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);

  if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Matrix"))) {
    args[1] = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@nrows"));
    args[0] = mrb_iv_get(mrb, other, mrb_intern_lit(mrb, "@ncols"));
    res = mrb_obj_new(mrb, mrb_class_get(mrb, "Matrix"), 2, args);
    mrb_matrix_get_data(mrb, other, &p_mat_other);
    mrb_matrix_get_data(mrb, res, &p_mat_res);

    if (p_mat->size2 != p_mat_other->size1) {
      mrb_raise(mrb, E_MATRIX_ERROR, "matrix dimensions don't match!");
    }
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, p_mat, p_mat_other, 0.0,
                   p_mat_res);
  } else if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Vector"))) {
    args[0] = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@nrows"));
    res = mrb_obj_new(mrb, mrb_class_get(mrb, "Vector"), 1, args);
    mrb_vector_get_data(mrb, other, &p_vec_other);
    mrb_vector_get_data(mrb, res, &p_vec_res);

    if (p_mat->size2 != p_vec_other->size) {
      mrb_raise(mrb, E_MATRIX_ERROR, "matrix dimensions don't match!");
    }
    gsl_blas_dgemv(CblasNoTrans, 1.0, p_mat, p_vec_other, 0.0, p_vec_res);
  }
  return res;
}

static mrb_value mrb_matrix_div(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_matrix *p_mat, *p_mat_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  mrb_matrix_get_data(mrb, other, &p_mat_other);
  if (p_mat->size1 != p_mat_other->size1 ||
      p_mat->size2 != p_mat_other->size2) {
    mrb_raise(mrb, E_MATRIX_ERROR, "matrix dimensions don't match!");
  }
  gsl_matrix_div_elements(p_mat, p_mat_other);
  return self;
}

static mrb_value mrb_matrix_transpose_self(mrb_state *mrb, mrb_value self) {
  gsl_matrix *p_mat;
  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (p_mat->size1 != p_mat->size2) {
    mrb_raise(mrb, E_MATRIX_ERROR, "matrix must be square!");
  }
  if (gsl_matrix_transpose(p_mat)) {
    mrb_raise(mrb, E_MATRIX_ERROR, "Cannot calculate transposed matrix");
  }
  return self;
}

static mrb_value mrb_matrix_transpose(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_matrix *p_mat, *p_mat_other;
  mrb_value args[2];
  // swap dimensions!
  args[1] = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@nrows"));
  args[0] = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@ncols"));
  other = mrb_obj_new(mrb, mrb_class_get(mrb, "Matrix"), 2, args);
  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  mrb_matrix_get_data(mrb, other, &p_mat_other);
  if (gsl_matrix_transpose_memcpy(p_mat_other, p_mat)) {
    mrb_raise(mrb, E_MATRIX_ERROR, "Cannot calculate transposed matrix");
  }
  return other;
}

static mrb_value mrb_matrix_swap_rows(mrb_state *mrb, mrb_value self) {
  gsl_matrix *p_mat;
  mrb_int i, j;
  mrb_get_args(mrb, "ii", &i, &j);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (gsl_matrix_swap_rows(p_mat, i, j)) {
    mrb_raise(mrb, E_MATRIX_ERROR, "Cannot swap rows");
  }
  return self;
}

static mrb_value mrb_matrix_swap_cols(mrb_state *mrb, mrb_value self) {
  gsl_matrix *p_mat;
  mrb_int i, j;
  mrb_get_args(mrb, "ii", &i, &j);

  // call utility for unwrapping @data into p_data:
  mrb_matrix_get_data(mrb, self, &p_mat);
  if (gsl_matrix_swap_columns(p_mat, i, j)) {
    mrb_raise(mrb, E_MATRIX_ERROR, "Cannot swap cols");
  }
  return self;
}

#pragma mark -
#pragma mark • Gem setup

void mrb_gsl_matrix_init(mrb_state *mrb) {
  struct RClass *gsl;

  mrb_load_string(mrb, "class MatrixError < Exception; end");

  gsl = mrb_define_class(mrb, "Matrix", mrb->object_class);
  mrb_define_method(mrb, gsl, "initialize", mrb_matrix_initialize,
                    MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "dup", mrb_matrix_dup, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "all", mrb_matrix_all, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "zero", mrb_matrix_zero, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "identity", mrb_matrix_identity, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "rnd_fill", mrb_matrix_rnd_fill, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "===", mrb_matrix_equal, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "[]", mrb_matrix_get_ij, MRB_ARGS_OPT(2));
  mrb_define_method(mrb, gsl, "row", mrb_matrix_get_row, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "col", mrb_matrix_get_col, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "[]=", mrb_matrix_set_ij, MRB_ARGS_REQ(3));
  mrb_define_method(mrb, gsl, "get_row", mrb_matrix_get_row, MRB_ARGS_OPT(1));
  mrb_define_method(mrb, gsl, "get_col", mrb_matrix_get_col, MRB_ARGS_OPT(1));
  mrb_define_method(mrb, gsl, "set_row", mrb_matrix_set_row, MRB_ARGS_REQ(2));
  mrb_define_method(mrb, gsl, "set_col", mrb_matrix_set_col, MRB_ARGS_REQ(2));

  mrb_define_method(mrb, gsl, "max", mrb_matrix_max, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "max_index", mrb_matrix_max_index,
                    MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "min", mrb_matrix_min, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "min_index", mrb_matrix_min_index,
                    MRB_ARGS_NONE());

  mrb_define_method(mrb, gsl, "add!", mrb_matrix_add, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "sub!", mrb_matrix_sub, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "mul!", mrb_matrix_mul, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "div!", mrb_matrix_div, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "^", mrb_matrix_prod, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "t!", mrb_matrix_transpose_self, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "t", mrb_matrix_transpose, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "swap_rows", mrb_matrix_swap_rows,
                    MRB_ARGS_REQ(2));
  mrb_define_method(mrb, gsl, "swap_cols", mrb_matrix_swap_cols,
                    MRB_ARGS_REQ(2));
}

max_stars_repo_name: UniTN-Mechatronics/mruby-gsl

max_issues_repo_name: UniTN-Mechatronics/mruby-gsl

max_forks_repo_name: UniTN-Mechatronics/mruby-gsl


text:
//
// gemmtest.c
//
// J. Makino
//    Time-stamp: <10/10/17 13:43:33 makino>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <timerlib.h>
#ifndef NOBLAS
#include <cblas.h>
#endif
#ifdef USEGDR
#include "gdrdgemm.h"
#endif

#define FTYPE double
#include <emmintrin.h>
typedef double v2df __attribute__((vector_size(16)));
typedef union {v2df v; double s[2];}v2u;

#ifndef USEGDR
void  gdrsetboardid(int boardid)
{}
#endif

void matmul2_host(int n,
		  FTYPE a[n][n],
		  FTYPE b[n][n],
		  FTYPE c[n][n])

{
    int i, j, k;
    for(i=0;i<n;i++){
	for(j=0;j<n;j++){
	    c[i][j]=0.0e0;
	    for(k=0;k<n;k++) c[i][j]+= a[i][k]*b[k][j];
	}
    }
}

// simplest version
void matmul_for_small_nk_0(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{
    // simplest version
    int i,j,k;
    for(j=0;j<n;j++)
	for(i=0;i<m;i++)
	    for(k=0;k<kk;k++)
		c[i][j] -= a[i][k]*b[k][j];
}

// XMM registers
#define X0	"%xmm0"
#define X1	"%xmm1"
#define X2	"%xmm2"
#define X3	"%xmm3"
#define X4	"%xmm4"
#define X5	"%xmm5"
#define X6	"%xmm6"
#define X7	"%xmm7"
#define X8		"%xmm8"
#define X9		"%xmm9"
#define X10		"%xmm10"
#define X11		"%xmm11"
#define X12		"%xmm12"
#define X13		"%xmm13"
#define X14		"%xmm14"
#define X15		"%xmm15"

#define LOADPD(mem, reg) asm("movapd %0, %"reg::"m"(mem));
#define STORPD(reg, mem) asm("movapd %"reg " , %0"::"m"(mem));
#define MOVNTPD(reg, mem) asm("movntpd %"reg " , %0"::"m"(mem));
#define MOVAPD(src, dst) asm("movapd " src "," dst);
#define MOVQ(src, dst) asm("movq " src "," dst);
#define BCAST0(reg) asm("shufpd $0x00, " reg ","  reg);
#define BCAST1(reg) asm("shufpd $0xff, " reg ","  reg);
#define MULPD(src, dst) asm("mulpd " src "," dst);
#define ADDPD(src, dst) asm("addpd " src ","  dst);
#define SUBPD(src, dst) asm("subpd "  src "," dst);

void daxpy(v2df a[], v2df b[], v2df c[], int n)
{
    int i;
    for(i=0;i<n; i+=8){
	b[i] += a[i]*c[0];
	b[i+1] += a[i+1]*c[0];
	b[i+2] += a[i+2]*c[0];
	b[i+3] += a[i+3]*c[0];
	b[i+4] += a[i+4]*c[0];
	b[i+5] += a[i+5]*c[0];
	b[i+6] += a[i+6]*c[0];
	b[i+7] += a[i+7]*c[0];
    }
}

void matmul_for_nk16_0a(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 16;
    int nh = n/2;
    register int k;
    v2df bcopy2[nh][kk];
    for(j=0;j<nh;j++)
	for(k=0;k<kk;k++)
	    //		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
	    bcopy2[j][k] = *(((v2df*)b)+j*n2+k);
    for(i=0;i<m;i+=2){
	//	BEGIN_TSC;
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));
	v2df *app = (v2df*) a[i+1];
	v2df * cpp = (v2df*) (&(c[i+1][0]));

	int k;

	for(j=0;j<nh;j+=2){
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cpp[j],X14);
	    LOADPD(cpp[j+1],X15);

	    for(k=0;k<8;k++){

		LOADPD(ap[k],X0);
		LOADPD(app[k],X2);
		LOADPD(bvp0[k*2],X4);
		LOADPD(bvp1[k*2],X5);
		LOADPD(bvp0[k*2+1],X8);
		LOADPD(bvp1[k*2+1],X9);
		MOVAPD(X4,X6);
		MOVAPD(X5,X7);
		MOVAPD(X8,X10);
		MOVAPD(X9,X11);
		MOVAPD(X0,X1);
		BCAST0(X0);
		BCAST1(X1);
		MOVAPD(X2,X3);
		BCAST0(X2);
		BCAST1(X3);
		MULPD(X0,X4);
		MULPD(X0,X5);
		MULPD(X1,X8);
		MULPD(X1,X9);
		SUBPD(X4,X12);
		SUBPD(X5,X13);
		SUBPD(X8,X12);
		SUBPD(X9,X13);
		MULPD(X2,X6);
		MULPD(X2,X7);
		MULPD(X3,X10);
		MULPD(X3,X11);
		SUBPD(X6,X14);
		SUBPD(X7,X15);
		SUBPD(X10,X14);
		SUBPD(X11,X15);
	    }

	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cpp[j+0]);
	    STORPD(X15,cpp[j+1]);
	}
    }
}
void matmul_for_nk16_test1(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j,k;
    int kk = 16;
    int nh = n/2;
    v2df bcopy[nh][kk];
    for(j=0;j<nh;j++)
	for(k=0;k<kk;k++)
		bcopy[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i+=2){

	v2df csub[nh][2];
	for(j=0;j<2;j++){
	    v2df * src = (v2df*) c[i+j];
	    for(k=0;k<nh;k++){
		csub[k][j]=src[k];
	    }
	}

	for(k=0;k<kk;k+=4){
	    v2df *ap = (v2df*) a[i]+k;
	    v2df *app = (v2df*) a[i+1]+k;
	    LOADPD(ap[0],X8);
	    LOADPD(ap[1],X9);
	    LOADPD(app[0],X10);
	    LOADPD(app[1],X11);
	    MOVAPD(X8,X0);
	    BCAST0(X0);
	    MOVAPD(X8,X1);
	    BCAST1(X1);
	    MOVAPD(X9,X2);
	    BCAST0(X2);
	    MOVAPD(X9,X3);
	    BCAST1(X3);
	    MOVAPD(X10,X4);
	    BCAST0(X4);
	    MOVAPD(X10,X5);
	    BCAST1(X5);
	    MOVAPD(X11,X6);
	    BCAST0(X6);
	    MOVAPD(X11,X7);
	    BCAST1(X7);
	    // 2x4 a, size doubled to 4x4
	    for(j=0;j<nh;j++){
		v2df * bvp = bcopy[j]+k;
		v2df * cvp = csub[k];
		LOADPD(bvp[0],X8);
		LOADPD(bvp[1],X9);

		MOVAPD(X8,X14);
		MOVAPD(X9,X15);
		MULPD(X0,X14);
		MULPD(X1,X15);

		LOADPD(bvp[2],X10);
		LOADPD(bvp[3],X11);

		LOADPD(cvp[0],X12);
		LOADPD(cvp[1],X13);

		SUBPD(X14,X12);
		SUBPD(X15,X12);
		MOVAPD(X10,X14);
		MOVAPD(X11,X15);
		MULPD(X2,X14);
		MULPD(X3,X15);
		SUBPD(X14,X12);
		SUBPD(X15,X12);
		MULPD(X4,X8);
		MULPD(X5,X9);
		MULPD(X6,X10);
		MULPD(X7,X11);
		SUBPD(X8,X13);
		SUBPD(X9,X13);
		SUBPD(X10,X13);
		SUBPD(X11,X13);
		STORPD(X12, cvp[0]);
		STORPD(X13, cvp[1]);
	    }
	}
	// copyback cmat
	for(j=0;j<2;j++){
	    v2df * dest = (v2df*) c[i+j];
	    for(k=0;k<nh;k++){
		dest[k] = csub[k][j];
	    }
	}
    }
}

void matmul_for_nk16_0c(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 16;
    int nh = n/2;
    register int k;
    v2df bcopy2[nh][kk];
#ifdef     PREFETCHL
#undef PREFETCHL
#endif
#define PREFETCHL 16

    for(i=0;i<PREFETCHL;i++){
	__builtin_prefetch((double*)a[i],0,0);
	__builtin_prefetch((double*)a[i]+8,0,0);
	__builtin_prefetch(c[i+8],1,0);
	__builtin_prefetch(c[i+8]+8,1,0);
    }
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));

	__builtin_prefetch((double*)a[i+PREFETCHL],0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+8,0,0);

	int k;

	for(j=0;j<nh;j+=4){
	    __builtin_prefetch(c[i+PREFETCHL]+j,1,0);
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    v2df * bvp2 = bcopy2[j+2];
	    v2df * bvp3 = bcopy2[j+3];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cp[j+2],X14);
	    LOADPD(cp[j+3],X15);

	    for(k=0;k<8;k++){

		LOADPD(ap[k],X0);
		LOADPD(bvp0[k*2],X4);
		LOADPD(bvp1[k*2],X5);
		LOADPD(bvp2[k*2],X6);
		LOADPD(bvp3[k*2],X7);
		MOVAPD(X0,X1);
		BCAST0(X0);
		BCAST1(X1);
		LOADPD(bvp0[k*2+1],X8);
		LOADPD(bvp1[k*2+1],X9);
		LOADPD(bvp2[k*2+1],X10);
		LOADPD(bvp3[k*2+1],X11);
		MULPD(X0,X4);
		MULPD(X0,X5);
		MULPD(X0,X6);
		MULPD(X0,X7);
		MULPD(X1,X8);
		MULPD(X1,X9);
		MULPD(X1,X10);
		MULPD(X1,X11);
		SUBPD(X4,X12);
		SUBPD(X5,X13);
		SUBPD(X6,X14);
		SUBPD(X7,X15);
		SUBPD(X8,X12);
		SUBPD(X9,X13);
		SUBPD(X10,X14);
		SUBPD(X11,X15);
	    }

	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cp[j+2]);
	    STORPD(X15,cp[j+3]);
	}
    }
}
void matmul_for_nk32_0(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 32;
    int nh = n/2;
    register int k;
    v2df bcopy2[nh][kk];
#ifdef     PREFETCHL
#undef PREFETCHL
#endif
#define PREFETCHL 8

    for(i=0;i<PREFETCHL;i++){
	__builtin_prefetch((double*)a[i],0,0);
	__builtin_prefetch((double*)a[i]+8,0,0);
	__builtin_prefetch((double*)a[i]+16,0,0);
	__builtin_prefetch((double*)a[i]+24,0,0);
	__builtin_prefetch(c[i+8],1,0);
	__builtin_prefetch(c[i+8]+8,1,0);
	__builtin_prefetch(c[i+8]+16,1,0);
	__builtin_prefetch(c[i+8]+24,1,0);
    }
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));

	__builtin_prefetch((double*)a[i+PREFETCHL],0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+8,0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+16,0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+24,0,0);

	int k;

	for(j=0;j<nh;j+=4){
	    __builtin_prefetch(c[i+PREFETCHL]+j,1,0);
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    v2df * bvp2 = bcopy2[j+2];
	    v2df * bvp3 = bcopy2[j+3];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cp[j+2],X14);
	    LOADPD(cp[j+3],X15);

	    for(k=0;k<16;k++){

		LOADPD(ap[k],X0);
		LOADPD(bvp0[k*2],X4);
		LOADPD(bvp1[k*2],X5);
		LOADPD(bvp2[k*2],X6);
		LOADPD(bvp3[k*2],X7);
		MOVAPD(X0,X1);
		BCAST0(X0);
		BCAST1(X1);
		MULPD(X0,X4);
		MULPD(X0,X5);
		MULPD(X0,X6);
		MULPD(X0,X7);
		LOADPD(bvp0[k*2+1],X8);
		LOADPD(bvp1[k*2+1],X9);
		LOADPD(bvp2[k*2+1],X10);
		LOADPD(bvp3[k*2+1],X11);
		MULPD(X1,X8);
		MULPD(X1,X9);
		MULPD(X1,X10);
		MULPD(X1,X11);
		SUBPD(X4,X12);
		SUBPD(X5,X13);
		SUBPD(X6,X14);
		SUBPD(X7,X15);
		SUBPD(X8,X12);
		SUBPD(X9,X13);
		SUBPD(X10,X14);
		SUBPD(X11,X15);
	    }

	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cp[j+2]);
	    STORPD(X15,cp[j+3]);
	}
    }
}

void matmul_for_nk16_0b(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 16;
    int nh = n/2;
    register int k;
    v2df bcopy2[nh][kk];
#ifdef     PREFETCHL
#undef PREFETCHL
#endif
#define PREFETCHL 16

    for(i=0;i<PREFETCHL;i++){
	__builtin_prefetch((double*)a[i],0,0);
	__builtin_prefetch((double*)a[i]+8,0,0);
	__builtin_prefetch(c[i+8],1,0);
	__builtin_prefetch(c[i+8]+8,1,0);
    }
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));

	__builtin_prefetch((double*)a[i+PREFETCHL],0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+8,0,0);

	int k;

	for(j=0;j<nh;j+=4){
	    __builtin_prefetch(c[i+PREFETCHL]+j,1,0);
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    v2df * bvp2 = bcopy2[j+2];
	    v2df * bvp3 = bcopy2[j+3];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cp[j+2],X14);
	    LOADPD(cp[j+3],X15);

	    LOADPD(ap[0],X0);
	    LOADPD(bvp0[0],X4);
	    LOADPD(bvp1[0],X5);
	    LOADPD(bvp2[0],X6);
	    LOADPD(bvp3[0],X7);
	    LOADPD(bvp0[1],X8);
	    LOADPD(bvp1[1],X9);
	    LOADPD(bvp2[1],X10);
	    LOADPD(bvp3[1],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[1],X0);

	    LOADPD(bvp0[2],X4);
	    LOADPD(bvp1[2],X5);
	    LOADPD(bvp2[2],X6);
	    LOADPD(bvp3[2],X7);
	    LOADPD(bvp0[3],X8);
	    LOADPD(bvp1[3],X9);
	    LOADPD(bvp2[3],X10);
	    LOADPD(bvp3[3],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[2],X0);
	    LOADPD(bvp0[4],X4);
	    LOADPD(bvp1[4],X5);
	    LOADPD(bvp2[4],X6);
	    LOADPD(bvp3[4],X7);
	    LOADPD(bvp0[5],X8);
	    LOADPD(bvp1[5],X9);
	    LOADPD(bvp2[5],X10);
	    LOADPD(bvp3[5],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[3],X0);
	    LOADPD(bvp0[6],X4);
	    LOADPD(bvp1[6],X5);
	    LOADPD(bvp2[6],X6);
	    LOADPD(bvp3[6],X7);
	    LOADPD(bvp0[7],X8);
	    LOADPD(bvp1[7],X9);
	    LOADPD(bvp2[7],X10);
	    LOADPD(bvp3[7],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[4],X0);

	    LOADPD(bvp0[8],X4);
	    LOADPD(bvp1[8],X5);
	    LOADPD(bvp2[8],X6);
	    LOADPD(bvp3[8],X7);
	    LOADPD(bvp0[9],X8);
	    LOADPD(bvp1[9],X9);
	    LOADPD(bvp2[9],X10);
	    LOADPD(bvp3[9],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);
	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[5],X0);
	    LOADPD(bvp0[10],X4);
	    LOADPD(bvp1[10],X5);
	    LOADPD(bvp2[10],X6);
	    LOADPD(bvp3[10],X7);
	    LOADPD(bvp0[11],X8);
	    LOADPD(bvp1[11],X9);
	    LOADPD(bvp2[11],X10);
	    LOADPD(bvp3[11],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[6],X0);
	    LOADPD(bvp0[12],X4);
	    LOADPD(bvp1[12],X5);
	    LOADPD(bvp2[12],X6);
	    LOADPD(bvp3[12],X7);
	    LOADPD(bvp0[13],X8);
	    LOADPD(bvp1[13],X9);
	    LOADPD(bvp2[13],X10);
	    LOADPD(bvp3[13],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[7],X0);
	    LOADPD(bvp0[14],X4);
	    LOADPD(bvp1[14],X5);
	    LOADPD(bvp2[14],X6);
	    LOADPD(bvp3[14],X7);
	    LOADPD(bvp0[15],X8);
	    LOADPD(bvp1[15],X9);
	    LOADPD(bvp2[15],X10);
	    LOADPD(bvp3[15],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cp[j+2]);
	    STORPD(X15,cp[j+3]);
	}
    }
}

void matmul_for_nk8_0d(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df awork[4];
    v2df awork2[4];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	double *ap=a[i];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));

	v2df * aa = (v2df*)(ap);
	__builtin_prefetch((double*)a[i+8],0,0);

	v2df acopy0=(v2df){a[i][0], a[i][0]};
	v2df acopy1=(v2df){a[i][1], a[i][1]};
	v2df acopy2=(v2df){a[i][2], a[i][2]};
	v2df acopy3=(v2df){a[i][3], a[i][3]};
	v2df acopy4=(v2df){a[i][4], a[i][4]};
	v2df acopy5=(v2df){a[i][5], a[i][5]};
	v2df acopy6=(v2df){a[i][6], a[i][6]};
	v2df acopy7=(v2df){a[i][7], a[i][7]};
	v2df zero=(v2df){0.0, 0.0};

	LOADPD(acopy0,X0);
	LOADPD(acopy1,X1);
	LOADPD(acopy2,X2);
	LOADPD(acopy3,X3);
	LOADPD(acopy4,X4);
	LOADPD(acopy5,X5);
	LOADPD(acopy6,X6);
	LOADPD(acopy7,X7);

	for(j=0;j<nh;j++){
	    __builtin_prefetch(c[i+8]+j,1,0);
	    v2df * bvp = bcopy2[j];
	    LOADPD(cp[j],X14);
	    LOADPD(bvp[0],X8);
	    LOADPD(bvp[1],X9);
	    MULPD(X0,X8);
	    MULPD(X1,X9);
	    LOADPD(bvp[2],X10);
	    LOADPD(bvp[3],X11);
	    ADDPD(X9,X8);
	    MULPD(X2,X10);
	    MULPD(X3,X11);
	    ADDPD(X11,X10);
	    LOADPD(bvp[4],X9);
	    LOADPD(bvp[5],X11);
	    LOADPD(bvp[6],X12);
	    LOADPD(bvp[7],X13);
	    MULPD(X4,X9);
	    MULPD(X5,X11);
	    ADDPD(X10,X8);
	    ADDPD(X11,X9);
	    MULPD(X6,X12);
	    MULPD(X7,X13);
	    ADDPD(X13,X12);
	    ADDPD(X9,X8);
	    ADDPD(X12,X8);
	    SUBPD(X8,X14);
	    STORPD(X14,cp[j]);
	}
    }
}

void matmul_for_nk8_0c(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df awork[4];
    v2df awork2[4];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	double *ap=a[i];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));

	v2df * aa = (v2df*)(ap);
	__builtin_prefetch((double*)a[i+8],0,0);

	register v2df acopy0=(v2df){a[i][0], a[i][0]};
	register v2df acopy1=(v2df){a[i][1], a[i][1]};
	register v2df acopy2=(v2df){a[i][2], a[i][2]};
	register v2df acopy3=(v2df){a[i][3], a[i][3]};
	register v2df acopy4=(v2df){a[i][4], a[i][4]};
	register v2df acopy5=(v2df){a[i][5], a[i][5]};
	register v2df acopy6=(v2df){a[i][6], a[i][6]};
	register v2df acopy7=(v2df){a[i][7], a[i][7]};

	for(j=0;j<nh;j++){
	    tmp =  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];

	    v2df * bp = bcopy2[j];
	    __builtin_prefetch(c[i+4]+j,1,0);
	    v2df *bvp = bp;
	    tmp += acopy0*bvp[0];
	    tmp +=acopy1*bvp[1];
	    tmp +=acopy2*bvp[2];
	    tmp +=acopy3*bvp[3];
	    tmp +=acopy4*bvp[4];
	    tmp +=acopy5*bvp[5];
	    tmp +=acopy6*bvp[6];
	    tmp +=acopy7*bvp[7];
	    cp[j] = ctmp -tmp;
	}
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);

}
void matmul_for_nk8_0b(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df awork[4];
    v2df awork2[4];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	double *ap=a[i];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));

	v2df * aa = (v2df*)(ap);
	__builtin_prefetch((double*)a[i+8],0,0);

	acopy[0]=(v2df){a[i][0], a[i][0]};
	acopy[1]=(v2df){a[i][1], a[i][1]};
	acopy[2]=(v2df){a[i][2], a[i][2]};
	acopy[3]=(v2df){a[i][3], a[i][3]};
	acopy[4]=(v2df){a[i][4], a[i][4]};
	acopy[5]=(v2df){a[i][5], a[i][5]};
	acopy[6]=(v2df){a[i][6], a[i][6]};
	acopy[7]=(v2df){a[i][7], a[i][7]};

	for(j=0;j<nh;j++){
	    tmp = tmp2=  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];

	    v2df * bp = bcopy2[j];
	    __builtin_prefetch(c[i+4]+j,1,0);
	    v2df *avp = acopy;
	    v2df *bvp = bp;
	    tmp += avp[0]*bvp[0];
	    tmp +=avp[1]*bvp[1];
	    tmp +=avp[2]*bvp[2];
	    tmp +=avp[3]*bvp[3];
	    tmp += avp[4]*bvp[4];
	    tmp +=avp[5]*bvp[5];
	    tmp +=avp[6]*bvp[6];
	    tmp +=avp[7]*bvp[7];
	    cp[j] = ctmp -tmp;
	}
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);

}
void matmul_for_nk8_0a(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df awork[4];
    v2df awork2[4];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i+=2){
	//	BEGIN_TSC;
	double *ap=a[i];
	double *ap2=a[i+1];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));
	v2df * cp2 = (v2df*) (&(c[i+1][0]));

	v2df * aa = (v2df*)(ap);
	__builtin_prefetch((double*)a[i+8],0,0);
	__builtin_prefetch((double*)a[i+9],0,0);

	acopy[0]=(v2df){a[i][0], a[i][0]};
	acopy[1]=(v2df){a[i][1], a[i][1]};
	acopy[2]=(v2df){a[i][2], a[i][2]};
	acopy[3]=(v2df){a[i][3], a[i][3]};
	acopy[4]=(v2df){a[i][4], a[i][4]};
	acopy[5]=(v2df){a[i][5], a[i][5]};
	acopy[6]=(v2df){a[i][6], a[i][6]};
	acopy[7]=(v2df){a[i][7], a[i][7]};

	aa = (v2df*)(ap2);
	acopy2[0]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	acopy2[1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	aa++;
	acopy2[2]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	acopy2[3]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	aa++;
	acopy2[4]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	acopy2[5]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	aa++;
	acopy2[6]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	acopy2[7]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	aa++;
	for(j=0;j<nh;j++){
	    tmp = tmp2=  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];
	    v2df ctmp2 = cp2[j] ;

	    v2df * bp = bcopy2[j];
	    __builtin_prefetch(c[i+4]+j,1,0);
	    __builtin_prefetch(c[i+5]+j,1,0);
	    v2df *avp = acopy;
	    v2df *avp2 = acopy2;
	    v2df *bvp = bp;
	    tmp += avp[0]*bvp[0];
	    tmp2 += avp2[0]*bvp[0];
	    tmp +=avp[1]*bvp[1];
	    tmp2+=avp2[1]*bvp[1];
	    tmp +=avp[2]*bvp[2];
	    tmp2+=avp2[2]*bvp[2];
	    tmp +=avp[3]*bvp[3];
	    tmp2+=avp2[3]*bvp[3];
	    tmp += avp[4]*bvp[4];
	    tmp2 += avp2[4]*bvp[4];
	    tmp +=avp[5]*bvp[5];
	    tmp2+=avp2[5]*bvp[5];
	    tmp +=avp[6]*bvp[6];
	    tmp2+=avp2[6]*bvp[6];
	    tmp +=avp[7]*bvp[7];
	    tmp2+=avp2[7]*bvp[7];
	    cp[j] = ctmp -tmp;
	    cp2[j] = ctmp2 -tmp2;

	}

    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);

}
void matmul_for_nk8_1(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j,k;
    const int kk = 8;
    const int kh = kk/2;
    int nh = n/2;
    v2df bcopy[n][kh];
    v2df acopy[kk][kh];
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(j=0;j<n;j++)
	for(k=0;k<kh;k++)
	    bcopy[j][k] = (v2df){b[k*2][j],b[k*2+1][j]};
    //    printf("copy b end\n");
    for(i=0;i<m;i+=kk){
	for(k=0;k<kk;k++){
	    v2df  *ak = (v2df*)(a[i+k]);
	    v2df * awp =acopy+k;
	    awp[0]=ak[0];
	    awp[1]=ak[1];
	    awp[2]=ak[2];
	    awp[3]=ak[3];
	}
	//	printf("copy a end\n");
	for(k=0;k<kk;k++){
	    v2u tmp, tmp1;
	    v2df * ap = acopy[k];
	    for(j=0;j<n;j+=2){
		tmp.v = ap[0]*bcopy[j][0]
		    + ap[1]*bcopy[j][1]
		    + ap[2]*bcopy[j][2]
		    + ap[3]*bcopy[j][3];
		tmp1.v = ap[0]*bcopy[j+1][0]
		    + ap[1]*bcopy[j+1][1]
		    + ap[2]*bcopy[j+1][2]
		    + ap[3]*bcopy[j+1][3];
		c[k+i][j] -= tmp.s[0]+tmp.s[1];
		c[k+i][j+1] -= tmp1.s[0]+tmp1.s[1];
	    }
	}
	//	printf("calc c end\n");
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);
}
void matmul_for_nk8_2(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j,k;
    const int kk = 8;
    const int kh = kk/2;
    int nh = n/2;
    v2df bcopy[nh][kk];
    v2df acopy[kk][kh];
    v2df ccopy[kk][kh];
    v2df acopy2[kk][kk];
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy[j][k] = *((v2df*)(b[k]+j+j));
    //    printf("copy b end\n");
    for(i=0;i<m;i+=kk){
	for(k=0;k<kk;k++){
	    __builtin_prefetch(a+i+k+8,0,0);
	    __builtin_prefetch(c+i+k+8,1,0);
	    v2df  *ak = (v2df*)(a[i+k]);
	    v2df * awp = (v2df*)(acopy+k);
	    v2df  *ck = (v2df*)(c[i+k]);
	    v2df * cwp = (v2df*)(ccopy+k);
	    awp[0]=ak[0];
	    awp[1]=ak[1];
	    awp[2]=ak[2];
	    awp[3]=ak[3];
	    cwp[0]=ck[0];
	    cwp[1]=ck[1];
	    cwp[2]=ck[2];
	    cwp[3]=ck[3];
	}
	for (j=0;j<n;j++){
	    double * ap = (double*)( acopy+j);
	    for (k=0;k<kk;k++){
		acopy2[j][k]=(v2df){ap[k],ap[k]};
	    }
	}
	//	printf("copy a end\n");
	for(k=0;k<kk;k++){
	    v2df * cp = (v2df*) ccopy[k];
	    v2df * ap = acopy2[k];
	    for(j=0;j<nh;j++){
		v2df * bp = bcopy[j];
		cp[j] -= ap[0]*bp[0]
		    + ap[1]*bp[1]
		    + ap[2]*bp[2]
		    + ap[3]*bp[3]
		    + ap[4]*bp[4]
		    + ap[5]*bp[5]
		    + ap[6]*bp[6]
		    + ap[7]*bp[7];
	    }
	}
	for(k=0;k<kk;k++){
	    v2df  *ck = (v2df*)(c[i+k]);
	    v2df * cwp = (v2df*)(ccopy+k);
#if 0
	    ck[0] = cwp[0];
	    ck[1] = cwp[1];
	    ck[2] = cwp[2];
	    ck[3] = cwp[3];
#endif
	    __builtin_ia32_movntpd((double*)(ck),cwp[0]);
	    __builtin_ia32_movntpd((double*)(ck+1),cwp[1]);
	    __builtin_ia32_movntpd((double*)(ck+2),cwp[2]);
	    __builtin_ia32_movntpd((double*)(ck+3),cwp[3]);
	}
	//	printf("calc c end\n");
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);
}
void matmul_for_nk8(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j,k;
    const int kk = 8;
    const int kh = kk/2;
    int nh = n/2;
    v2df bcopy[nh][kk];
    v2df acopy2[kk][kk];
    //    unsigned long bpcount, apcount, dotcount;
    //    bpcount= apcount= dotcount=0;

    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy[j][k] = *((v2df*)(b[k]+j+j));
    //    END_TSC(bpcount);
    //    printf("copy b end\n");
#pragma omp parallel for private(i,j,k,acopy2)	 schedule(static)
    for(i=0;i<m;i+=kk){
	//	BEGIN_TSC;
	for(k=0;k<kk;k++){
	    __builtin_prefetch(a+i+k+16,0,0);
	    __builtin_prefetch(c+i+k+16,1,0);
	}
	for (j=0;j<n;j++){
	    double * ap = (double*)( a[i+j]);
	    for (k=0;k<kk;k++){
		acopy2[j][k]=(v2df){ap[k],ap[k]};
	    }
	}
	//	END_TSC(apcount);
	//	printf("copy a end\n");
	//	BEGIN_TSC;
	for(k=0;k<kk;k++){
	    v2df * cp = (v2df*) (c[i+k]);
	    v2df * ap = acopy2[k];
	    for(j=0;j<nh;j++){
		v2df * bp = bcopy[j];
		cp[j] -= ap[0]*bp[0] + ap[1]*bp[1]
		    + ap[2]*bp[2]    + ap[3]*bp[3]
		    + ap[4]*bp[4]    + ap[5]*bp[5]
		    + ap[6]*bp[6]    + ap[7]*bp[7];
	    }
	}
	//	printf("calc c end\n");
	//	END_TSC(dotcount);
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //   (double)bpcount, (double)apcount, (double)dotcount);
}

void matmul_for_nk16(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{
    int ii;
    int dm = (m+63)/64;
    dm*= 16;
#pragma omp parallel for private(ii)	 schedule(static)
    for(ii=0;ii<4;ii++){
	int ifirst, iend;
	ifirst = ii*dm;
	iend = ifirst+dm;
	if (iend > m) iend = m;
	//	fprintf(stderr, "m, i, ifirst, iend = %d %d %d %d\n", m, ii, ifirst, iend);
	if (ifirst < m){
	    matmul_for_nk16_0a(n1, a[ifirst], n2, b, n3, c[ifirst], iend-ifirst, n);
	}
    }

}

void matmul_for_nk32(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{
    int ii;
    int dm = (m+127)/128;
    dm*= 32;
#pragma omp parallel for private(ii)	 schedule(static)
    for(ii=0;ii<4;ii++){
	int ifirst, iend;
	ifirst = ii*dm;
	iend = ifirst+dm;
	if (iend > m) iend = m;
	//	fprintf(stderr, "m, i, ifirst, iend = %d %d %d %d\n", m, ii, ifirst, iend);
	if (ifirst < m){
	    matmul_for_nk32_0(n1, a[ifirst], n2, b, n3, c[ifirst], iend-ifirst, n);
	}
    }

}

#ifndef USEGDR
void gdrsetforceswapab(){}
void gdrresetforceswapab(){}
void gdrsetskipsendjmat(){};
void gdrresetskipsendjmat(){}
void gdrsetnboards(){}
void set_matmul_msg_level(int level){}
#endif

#define N 1024
#define K 16
int main()
{
    double a[N][K];
    double c[N][K];
    double b[K][K];
    int i,j;
    for(i=0;i<N;i++)
	for(j=0;j<K;j++){
	    a[i][j] = i*j;
	   c[i][j] = 0;
	}
    for(i=0;i<K;i++)
	for(j=0;j<K;j++)
	    b[i][j] = i*j;

    int k;
    unsigned long int start, end;
    rdtscl(&start);
#define NT 5000
    for(i=0;i<NT;i++){
	matmul_for_nk16_test1(K, a, K, b, K, c, N, K);
    }
    rdtscl(&end);
    printf("cycles = %e; %e ops/clock\n",
	   (double)(end-start), (N*K*K*(NT+0.0))/((double)(end-start)));
}

ext: c


text:
//This does CELL (~soma) stage of LSTM (long short-term memory) model.
//This requires each neuron to have 4 input time series, Xc, Xi, Xf, Xo,
//where Xc is the usual (or "cellular") input and Xi, Xf, Xo the inputs for the input, forget, output gates.
//Xc, Xi, Xf, Xo are the output of separate linear IN stages (weights and baises).

//For dim=0, C[:,t] = tanh{Xc[:,t] + Uc*Y[:,t-1]} \n";
//           I[:,t] = sig{Xi[:,t] + Ui*Y[:,t-1]} \n";
//           F[:,t] = sig{Xf[:,t] + Uf*Y[:,t-1]} \n";
//           O[:,t] = sig{Xo[:,t] + Uo*Y[:,t-1]} \n";
//           H[:,t] = I[:,t].*C[:,t] + F[:,t].*H[:,t-1] \n";
//           Y[:,t] = O[:,t].*tanh{H[:,t]} \n";
//with sizes Xc, Xi, Xf, Xo: N x T \n";
//           Uc, Ui, Uf, Uo: N x N \n";
//           Y             : N x T \n";
//
//For dim=1, C[t,:] = tanh{Xc[t,:] + Y[t-1,:]*Uc} \n";
//           I[t,:] = sig{Xi[t,:] + Y[t-1,:]*Ui} \n";
//           F[t,:] = sig{Xf[t,:] + Y[t-1,:]*Uf} \n";
//           O[t,:] = sig{Xo[t,:] + Y[t-1,:]*Uo} \n";
//           H[t,:] = I[t,:].*C[t,:] + F[t,:].*H[t-1,:] \n";
//           Y[t,:] = O[t,:].*tanh{H[t,:]} \n";
//with sizes Xc, Xi, Xf, Xo: T x N \n";
//           Uc, Ui, Uf, Uo: N x N \n";
//           Y             : T x N \n";
//
//where sig is the logistic (sigmoid) nonlinearity = 1/(1+exp(-x)),
//I is the input gate, F is the forget gate, O is the output gate,
//C is the "cell input activation vector",
//H is an intermediate (hidden) vector (sometimes called the "cell state vector"),
//Uc, Ui, Uf, Uo are NxN matrices, and Y is the final output (sometimes called the "hidden state vector").

//Note that, the neurons of a layer are independent only if Uc, Ui, Uf, Uo are diagonal matrices.
//This is only really a CELL (~soma) stage in that case.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cblas.h>

#ifdef I
    #undef I
#endif

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int lstm4_s (float *Y, const float *Xc, const float *Xi, const float *Xf, const float *Xo, const float *Uc, const float *Ui, const float *Uf, const float *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int lstm4_d (double *Y, const double *Xc, const double *Xi, const double *Xf, const double *Xo, const double *Uc, const double *Ui, const double *Uf, const double *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int lstm4_inplace_s (float *Xc, const float *Xi, const float *Xf, const float *Xo, const float *Uc, const float *Ui, const float *Uf, const float *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int lstm4_inplace_d (double *Xc, const double *Xi, const double *Xf, const double *Xo, const double *Uc, const double *Ui, const double *Uf, const double *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int lstm4_s (float *Y, const float *Xc, const float *Xi, const float *Xf, const float *Xo, const float *Uc, const float *Ui, const float *Uf, const float *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    size_t nT, tN;

    float *C, *I, *F, *O, *H;
    if (!(C=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(I=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(F=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(O=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_s: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        //C[0] = tanhf(Xc[0]);
        //I[0] = 1.0f / (1.0f+expf(-Xi[0]));
        //F[0] = 1.0f / (1.0f+expf(-Xf[0]));
        //O[0] = 1.0f / (1.0f+expf(-Xo[0]));
        //H[0] = I[0] * C[0];
        //Y[0] = O[0] * tanhf(H[0]);
        H[0] = tanhf(Xc[0]) / (1.0f+expf(-Xi[0]));
        Y[0] = tanhf(H[0]) / (1.0f+expf(-Xo[0]));
        for (size_t t=1; t<T; ++t)
        {
            //C[0] = tanhf(Xc[t]+Uc[0]*Y[t-1]);
            //I[0] = 1.0f / (1.0f+expf(-Xi[t]-Ui[0]*Y[t-1]));
            //F[0] = 1.0f / (1.0f+expf(-Xf[t]-Uf[0]*Y[t-1]));
            //O[0] = 1.0f / (1.0f+expf(-Xo[t]-Uo[0]*Y[t-1]));
            //H[0] = I[0]*C[0] + F[0]*H[0];
            //Y[t] = O[0] * tanhf(H[0]);
            H[0] = tanhf(Xc[t]+Uc[0]*Y[t-1])/(1.0f+expf(-Xi[t]-Ui[0]*Y[t-1])) + H[0]/(1.0f+expf(-Xf[t]-Uf[0]*Y[t-1]));
            Y[t] = tanhf(H[0]) / (1.0f+expf(-Xo[t]-Uo[0]*Y[t-1]));
        }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                //C[n] = tanhf(Xc[n]);
                //I[n] = 1.0f / (1.0f+expf(-Xi[n]));
                //F[n] = 1.0f / (1.0f+expf(-Xf[n]));
                //O[n] = 1.0f / (1.0f+expf(-Xo[n]));
                //H[n] = I[n] * C[n];
                //Y[n] = O[n] * tanhf(H[n]);
                H[n] = tanhf(Xc[n]) / (1.0f+expf(-Xi[n]));
                Y[n] = tanhf(H[n]) / (1.0f+expf(-Xo[n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_scopy((int)N,&Xc[tN],1,C,1); cblas_scopy((int)N,&Xi[tN],1,I,1);
                cblas_scopy((int)N,&Xf[tN],1,F,1); cblas_scopy((int)N,&Xo[tN],1,O,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Y[tN-N],1,o,C,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Y[tN-N],1,o,I,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Y[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    //C[n] = tanhf(C[n]);
                    //I[n] = 1.0f / (1.0f+expf(-I[n]));
                    //F[n] = 1.0f / (1.0f+expf(-F[n]));
                    //O[n] = 1.0f / (1.0f+expf(-O[n]));
                    //H[n] = I[n]*C[n] + F[n]*H[n];
                    //Y[tN+n] = O[n] * tanhf(H[n]);
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Y[tN+n] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                //C[n] = tanhf(Xc[nT])
                //I[n] = 1.0f / (1.0f+expf(-Xi[nT]));
                //F[n] = 1.0f / (1.0f+expf(-Xf[nT]));
                //O[n] = 1.0f / (1.0f+expf(-Xo[nT]));
                //H[n] = I[n] * C[n];
                //Y[nT] = O[n] * tanhf(H[n]);
                H[n] = tanhf(Xc[nT]) / (1.0f+expf(-Xi[nT]));
                Y[nT] = tanhf(H[n]) / (1.0f+expf(-Xo[nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&Xc[t],(int)T,C,1); cblas_scopy((int)N,&Xi[t],(int)T,I,1);
                cblas_scopy((int)N,&Xf[t],(int)T,F,1); cblas_scopy((int)N,&Xo[t],(int)T,O,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Y[t-1],(int)T,o,C,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Y[t-1],(int)T,o,I,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Y[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    //C[n] = tanhf(C[n]);
                    //I[n] = 1.0f / (1.0f+expf(-I[n]));
                    //F[n] = 1.0f / (1.0f+expf(-F[n]));
                    //O[n] = 1.0f / (1.0f+expf(-O[n]));
                    //H[n] = I[n]*C[n] + F[n]*H[n];
                    //Y[t+n*T] = O[n] * tanhf(H[n]);
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Y[t+n*T] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanhf(Xc[nT]) / (1.0f+expf(-Xi[nT]));
                Y[nT] = tanhf(H[n]) / (1.0f+expf(-Xo[nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&Xc[t],(int)T,C,1); cblas_scopy((int)N,&Xi[t],(int)T,I,1);
                cblas_scopy((int)N,&Xf[t],(int)T,F,1); cblas_scopy((int)N,&Xo[t],(int)T,O,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Y[t-1],(int)T,o,C,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Y[t-1],(int)T,o,I,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Y[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Y[t+n*T] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanhf(Xc[n]) / (1.0f+expf(-Xi[n]));
                Y[n] = tanhf(H[n]) / (1.0f+expf(-Xo[n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_scopy((int)N,&Xc[tN],1,C,1); cblas_scopy((int)N,&Xi[tN],1,I,1);
                cblas_scopy((int)N,&Xf[tN],1,F,1); cblas_scopy((int)N,&Xo[tN],1,O,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Y[tN-N],1,o,C,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Y[tN-N],1,o,I,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Y[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Y[tN+n] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lstm4_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int lstm4_d (double *Y, const double *Xc, const double *Xi, const double *Xf, const double *Xo, const double *Uc, const double *Ui, const double *Uf, const double *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    size_t nT, tN;

    double *C, *I, *F, *O, *H;
    if (!(C=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(I=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(F=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(O=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_d: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        H[0] = tanh(Xc[0]) / (1.0+exp(-Xi[0]));
        Y[0] = tanh(H[0]) / (1.0+exp(-Xo[0]));
        for (size_t t=1; t<T; ++t)
        {
            H[0] = tanh(Xc[t]+Uc[0]*Y[t-1])/(1.0+exp(-Xi[t]-Ui[0]*Y[t-1])) + H[0]/(1.0+exp(-Xf[t]-Uf[0]*Y[t-1]));
            Y[t] = tanh(H[0]) / (1.0+exp(-Xo[t]-Uo[0]*Y[t-1]));
        }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanh(Xc[n]) / (1.0+exp(-Xi[n]));
                Y[n] = tanh(H[n]) / (1.0+exp(-Xo[n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_dcopy((int)N,&Xc[tN],1,C,1); cblas_dcopy((int)N,&Xi[tN],1,I,1);
                cblas_dcopy((int)N,&Xf[tN],1,F,1); cblas_dcopy((int)N,&Xo[tN],1,O,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Y[tN-N],1,o,C,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Y[tN-N],1,o,I,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Y[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Y[tN+n] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanh(Xc[nT]) / (1.0+exp(-Xi[nT]));
                Y[nT] = tanh(H[n]) / (1.0+exp(-Xo[nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&Xc[t],(int)T,C,1); cblas_dcopy((int)N,&Xi[t],(int)T,I,1);
                cblas_dcopy((int)N,&Xf[t],(int)T,F,1); cblas_dcopy((int)N,&Xo[t],(int)T,O,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Y[t-1],(int)T,o,C,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Y[t-1],(int)T,o,I,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Y[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Y[t+n*T] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanh(Xc[nT]) / (1.0+exp(-Xi[nT]));
                Y[nT] = tanh(H[n]) / (1.0+exp(-Xo[nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&Xc[t],(int)T,C,1); cblas_dcopy((int)N,&Xi[t],(int)T,I,1);
                cblas_dcopy((int)N,&Xf[t],(int)T,F,1); cblas_dcopy((int)N,&Xo[t],(int)T,O,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Y[t-1],(int)T,o,C,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Y[t-1],(int)T,o,I,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Y[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Y[t+n*T] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanh(Xc[n]) / (1.0+exp(-Xi[n]));
                Y[n] = tanh(H[n]) / (1.0+exp(-Xo[n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_dcopy((int)N,&Xc[tN],1,C,1); cblas_dcopy((int)N,&Xi[tN],1,I,1);
                cblas_dcopy((int)N,&Xf[tN],1,F,1); cblas_dcopy((int)N,&Xo[tN],1,O,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Y[tN-N],1,o,C,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Y[tN-N],1,o,I,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Y[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Y[tN+n] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lstm4_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int lstm4_inplace_s (float *Xc, const float *Xi, const float *Xf, const float *Xo, const float *Uc, const float *Ui, const float *Uf, const float *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    size_t nT, tN;

    float *C, *I, *F, *O, *H;
    if (!(C=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(I=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(F=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(O=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in lstm4_inplace_s: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        H[0] = tanhf(Xc[0]) / (1.0f+expf(-Xi[0]));
        Xc[0] = tanhf(H[0]) / (1.0f+expf(-Xo[0]));
        for (size_t t=1; t<T; ++t)
        {
            H[0] = tanhf(Xc[t]+Uc[0]*Xc[t-1])/(1.0f+expf(-Xi[t]-Ui[0]*Xc[t-1])) + H[0]/(1.0f+expf(-Xf[t]-Uf[0]*Xc[t-1]));
            Xc[t] = tanhf(H[0]) / (1.0f+expf(-Xo[t]-Uo[0]*Xc[t-1]));
        }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanhf(Xc[n]) / (1.0f+expf(-Xi[n]));
                Xc[n] = tanhf(H[n]) / (1.0f+expf(-Xo[n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_scopy((int)N,&Xc[tN],1,C,1); cblas_scopy((int)N,&Xi[tN],1,I,1);
                cblas_scopy((int)N,&Xf[tN],1,F,1); cblas_scopy((int)N,&Xo[tN],1,O,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Xc[tN-N],1,o,C,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Xc[tN-N],1,o,I,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Xc[tN-N],1,o,F,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Xc[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Xc[tN+n] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanhf(Xc[nT]) / (1.0f+expf(-Xi[nT]));
                Xc[nT] = tanhf(H[n]) / (1.0f+expf(-Xo[nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&Xc[t],(int)T,C,1); cblas_scopy((int)N,&Xi[t],(int)T,I,1);
                cblas_scopy((int)N,&Xf[t],(int)T,F,1); cblas_scopy((int)N,&Xo[t],(int)T,O,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Xc[t-1],(int)T,o,C,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Xc[t-1],(int)T,o,I,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Xc[t-1],(int)T,o,F,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Xc[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Xc[t+n*T] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanhf(Xc[nT]) / (1.0f+expf(-Xi[nT]));
                Xc[nT] = tanhf(H[n]) / (1.0f+expf(-Xo[nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&Xc[t],(int)T,C,1); cblas_scopy((int)N,&Xi[t],(int)T,I,1);
                cblas_scopy((int)N,&Xf[t],(int)T,F,1); cblas_scopy((int)N,&Xo[t],(int)T,O,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Xc[t-1],(int)T,o,C,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Xc[t-1],(int)T,o,I,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Xc[t-1],(int)T,o,F,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Xc[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Xc[t+n*T] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanhf(Xc[n]) / (1.0f+expf(-Xi[n]));
                Xc[n] = tanhf(H[n]) / (1.0f+expf(-Xo[n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_scopy((int)N,&Xc[tN],1,C,1); cblas_scopy((int)N,&Xi[tN],1,I,1);
                cblas_scopy((int)N,&Xf[tN],1,F,1); cblas_scopy((int)N,&Xo[tN],1,O,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Xc[tN-N],1,o,C,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Xc[tN-N],1,o,I,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Xc[tN-N],1,o,F,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Xc[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanhf(C[n])/(1.0f+expf(-I[n])) + H[n]/(1.0f+expf(-F[n]));
                    Xc[tN+n] = tanhf(H[n]) / (1.0f+expf(-O[n]));
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lstm4_inplace_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int lstm4_inplace_d (double *Xc, const double *Xi, const double *Xf, const double *Xo, const double *Uc, const double *Ui, const double *Uf, const double *Uo, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    size_t nT, tN;

    double *C, *I, *F, *O, *H;
    if (!(C=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(I=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(F=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(O=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in lstm4_inplace_d: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        H[0] = tanh(Xc[0]) / (1.0+exp(-Xi[0]));
        Xc[0] = tanh(H[0]) / (1.0+exp(-Xo[0]));
        for (size_t t=1; t<T; ++t)
        {
            H[0] = tanh(Xc[t]+Uc[0]*Xc[t-1])/(1.0+exp(-Xi[t]-Ui[0]*Xc[t-1])) + H[0]/(1.0+exp(-Xf[t]-Uf[0]*Xc[t-1]));
            Xc[t] = tanh(H[0]) / (1.0+exp(-Xo[t]-Uo[0]*Xc[t-1]));
        }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanh(Xc[n]) / (1.0+exp(-Xi[n]));
                Xc[n] = tanh(H[n]) / (1.0+exp(-Xo[n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_dcopy((int)N,&Xc[tN],1,C,1); cblas_dcopy((int)N,&Xi[tN],1,I,1);
                cblas_dcopy((int)N,&Xf[tN],1,F,1); cblas_dcopy((int)N,&Xo[tN],1,O,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Xc[tN-N],1,o,C,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Xc[tN-N],1,o,I,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Xc[tN-N],1,o,F,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Xc[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Xc[tN+n] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanh(Xc[nT]) / (1.0+exp(-Xi[nT]));
                Xc[nT] = tanh(H[n]) / (1.0+exp(-Xo[nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&Xc[t],(int)T,C,1); cblas_dcopy((int)N,&Xi[t],(int)T,I,1);
                cblas_dcopy((int)N,&Xf[t],(int)T,F,1); cblas_dcopy((int)N,&Xo[t],(int)T,O,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uc,(int)N,&Xc[t-1],(int)T,o,C,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ui,(int)N,&Xc[t-1],(int)T,o,I,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Xc[t-1],(int)T,o,F,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uo,(int)N,&Xc[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Xc[t+n*T] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                H[n] = tanh(Xc[nT]) / (1.0+exp(-Xi[nT]));
                Xc[nT] = tanh(H[n]) / (1.0+exp(-Xo[nT]));
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&Xc[t],(int)T,C,1); cblas_dcopy((int)N,&Xi[t],(int)T,I,1);
                cblas_dcopy((int)N,&Xf[t],(int)T,F,1); cblas_dcopy((int)N,&Xo[t],(int)T,O,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Xc[t-1],(int)T,o,C,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Xc[t-1],(int)T,o,I,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Xc[t-1],(int)T,o,F,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Xc[t-1],(int)T,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Xc[t+n*T] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                H[n] = tanh(Xc[n]) / (1.0+exp(-Xi[n]));
                Xc[n] = tanh(H[n]) / (1.0+exp(-Xo[n]));
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N;
                cblas_dcopy((int)N,&Xc[tN],1,C,1); cblas_dcopy((int)N,&Xi[tN],1,I,1);
                cblas_dcopy((int)N,&Xf[tN],1,F,1); cblas_dcopy((int)N,&Xo[tN],1,O,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uc,(int)N,&Xc[tN-N],1,o,C,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ui,(int)N,&Xc[tN-N],1,o,I,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Xc[tN-N],1,o,F,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uo,(int)N,&Xc[tN-N],1,o,O,1);
                for (size_t n=0u; n<N; ++n)
                {
                    H[n] = tanh(C[n])/(1.0+exp(-I[n])) + H[n]/(1.0+exp(-F[n]));
                    Xc[tN+n] = tanh(H[n]) / (1.0+exp(-O[n]));
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in lstm4_inplace_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

//#define DEBUG
#ifndef DEBUG
#define IN 64
#define OUT 32
#define NUM_ITS 1000000
#else
#define IN 4
#define OUT 2
#define NUM_ITS 1
#endif
#define INIT_VAL 0.4711

#ifdef USE_BLAS
#include <cblas.h>
#endif

double *
mtrx(int rows, int cols) {
    int row, col;
    double *d, *it;
    d = (double *) malloc(rows * cols * sizeof(double));
    it = d;
    for (row = 0; row < rows; row++) {
        for (col = 0; col < cols; col++) {
            // Not very random, but makes debugging easier.
            *it = INIT_VAL;
            it++;
        }
    }
    return d;
}

void
forward(const int rows, const int cols, double *W, double *x,
        double *a) {
    #ifndef USE_BLAS
    int row, col;
    double *W_it, *x_it, *a_it;

    W_it = W;
    a_it = a;
    for (row = 0; row < rows; row++) {
        *a_it = 0.0;
        x_it = x;
        for (col = 0; col < cols; col++) {
            *a_it += *W_it * *x_it;
            W_it++;
            x_it++;
        }
        *a_it = tanh(*a_it);
        a_it++;
    }
    #else /* USE_BLAS */
    int row;

    cblas_dgemv(CblasColMajor, CblasNoTrans, rows, cols, 1.0, W, rows, x, 1,
            0.0, a, 1);
    for (row = 0; row < rows; row++) {
        a[row] = tanh(a[row]);
    }
    #endif /* USE_BLAS */
}

void
backward(const int rows, const int cols, double *x, double *a, double *d,
        double *b) {
    #ifndef USE_BLAS
    int row, col;
    double ad;
    double *x_it, *a_it, *d_it, *b_it;

    a_it = a;
    d_it = d;
    b_it = b;
    for (row = 0; row < rows; row++) {
        ad = (1.0 - pow(*a_it, 2)) * *d_it;

        x_it = x;
        for (col = 0; col < cols; col++) {
            *b_it = ad * *x_it;
            b_it++;
            x_it++;
        }
        a_it++;
        d_it++;
    }
    #else /* USE_BLAS */
    int row;
    double *tmp;

    tmp = (double *) malloc(rows * sizeof(double));
    for (row = 0; row < rows; row++) {
        tmp[row] = (1.0 - pow(a[row], 2)) * d[row];
    }
    cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, rows, cols, 1, 1.0,
            tmp, 1, x, 1, 0.0, b, rows);

    free(tmp);
    #endif /* USE_BLAS */
}

int
main() {
    int i;
    double *W, *x, *d, *a, *b;
    struct timespec tic, toc;

    W = mtrx(OUT, IN);
    x = mtrx(IN, 1);
    d = mtrx(OUT, 1);
    a = malloc(OUT * sizeof(double));
    b = malloc(OUT * IN * sizeof(double));

    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &tic);
    for (i = 0; i < NUM_ITS; i++) {
        forward(OUT, IN, W, x, a);
        backward(OUT, IN, x, a, d, b);
    }
    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &toc);

    #ifdef DEBUG
    for (int j = 0; j < OUT; j++) {
        printf("%f ", a[j]);
    }
    printf("\n");

    double *b_it = b;
    for (int u = 0; u < OUT; u++) {
        for (int v = 0; v < IN; v++) {
            printf("%f ", *b_it);
            b_it++;
        }
        printf("\n");
    }
    #endif

    printf("%f\n", (toc.tv_nsec - tic.tv_nsec) / 100000000.0);

    free(b);
    free(a);
    free(d);
    free(x);
    free(W);

    return EXIT_SUCCESS;
}

max_stars_repo_name: ninjin/ppod

max_issues_repo_name: ninjin/ppod

max_forks_repo_name: ninjin/ppod


text:
#ifndef SCOL_H
#define SCOL_H
#include "Common.h"
// SP_IVP.h
// Spherical Collapse IVP
//#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <cvode/cvode.h>               /* prototypes for CVODE fcts., consts.  */
#include <nvector/nvector_serial.h>    /* access to serial N_Vector            */
#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix            */
#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver      */
#include <cvode/cvode_direct.h>        /* access to CVDls interface            */
#include <sundials/sundials_types.h>   /* defs. of realtype, sunindextype      */
#include <sundials/sundials_math.h>    /* contains the macros ABS, SUNSQR, EXP */

#define Ith(v,i)    NV_Ith_S(v,i-1)         /* Ith numbers components 1..NEQ */
#define IJth(A,i,j) SM_ELEMENT_D(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */

#define NEQ     2                /* number of equations  */
#define Y1      RCONST(0.0)      /* initial y components */
#define RTOL    RCONST(1.0e-5)   /* scalar relative tolerance            */
#define ATOL1   RCONST(1.0e-11)   /* vector absolute tolerance components */
#define ATOL2   RCONST(1.0e-11)
#define T00      RCONST(-10.41431317630211772495840705232694745063781738281250)      /* initial time           */
#define NOUT    1024               /* number of output times */
#define NOUT_DC    10               /* number of output times */

#define HALF    RCONST(0.5)      /* 0.5                    */
#define ONE     RCONST(1.0)      /* 1.0                    */
#define TWO     RCONST(2.0)
#define THREE   RCONST(3.0)       /* 3.0                    */
#define DELTA1  RCONST(3e-5)   /* Delta_i1 */
#define DELTA2  RCONST(0.00012)   /* Delta_i2 */
#define EPSILON RCONST(1.0e-9)
#define ZERO    RCONST(0.0)
#define NINE    RCONST(9.0)
#define TEN     RCONST(10.0)
#define Gnewton RCONST(4.302e-09)
#define coef    RCONST(8987404.41) // 1/H0^2

// structures to store spherical collapse calculations
typedef struct arrays{
  int count;
  double xx[1002];
  double yy[1002]; } *arrays_T;

typedef struct arrays3D{
  int count;
  double xx[1002];
  double yy[1002];
  double zz[1002]; } *arrays_T3;

typedef struct usdat {
  realtype IC;
  realtype OM;
  realtype Rth;
  realtype T1;
  double par1;
  double par2;
  double par3;
  int mymg;
  double maxt;
  gsl_spline *spline;
  gsl_interp_accel *acc;
} *UserData;

extern int check_flagscol(void *flagvalue, const char *funcname, int opt); //

class SCOL {
public:

// used functions in example file
   double maxP_zeta(double sig2, double dsig2dR, double OM, double Z);
   double Delta_Lambda(double OM, double Z); //
   // solves for y_enviornment
   int yenv(double OM_REAL, double XF, double delta_envi, arrays_T xxyy); // gives the environmental dependence of spherical collapse
   // solves for y_halo
   int SphericalCollapse(double *dC, arrays_T3 xxyyzz, UserData data_vec, double TMULT_REAL, double delta_g); // spherical collapse solver
   // solves for a_virial
   double myscol(double myscolparams[], double acol, double omega0, double Rthp, double sig1, double sig2, double pars[], int mymg, int yenvf); // solves for virial quantities and stores them in array myscolparams

   void PrintOutput(realtype t, realtype y1, realtype y2);
   void PrintRootInfo(int root_f1); //

   double funcscol(double xi, void *user_data); //

};

#endif

ext: h

max_stars_repo_name: PedroCarrilho/ReACT

max_issues_repo_name: PedroCarrilho/ReACT

max_forks_repo_name: PedroCarrilho/ReACT


text:
/* interpolation/accel.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>

gsl_interp_accel *
gsl_interp_accel_alloc (void)
{
  gsl_interp_accel *a = (gsl_interp_accel *) malloc (sizeof (gsl_interp_accel));
  if (a == 0)
    {
      GSL_ERROR_NULL("could not allocate space for gsl_interp_accel", GSL_ENOMEM);
    }

  a->cache = 0;
  a->hit_count = 0;
  a->miss_count = 0;

  return a;
}

int
gsl_interp_accel_reset (gsl_interp_accel * a)
{
  a->cache = 0;
  a->hit_count = 0;
  a->miss_count = 0;

  return GSL_SUCCESS;
}

#ifndef HIDE_INLINE_STATIC
size_t
gsl_interp_accel_find (gsl_interp_accel * a, const double xa[], size_t len, double x)
{
  size_t x_index = a->cache;

  if (x < xa[x_index])
    {
      a->miss_count++;
      a->cache = gsl_interp_bsearch (xa, x, 0, x_index);
    }
  else if (x > xa[x_index + 1])
    {
      a->miss_count++;
      a->cache = gsl_interp_bsearch (xa, x, x_index, len - 1);
    }
  else
    {
      a->hit_count++;
    }

  return a->cache;
}
#endif

void
gsl_interp_accel_free (gsl_interp_accel * a)
{
  free (a);
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include<stdio.h>
#include<sys/time.h>
#include <stdlib.h>

#include <cblas.h>

#include "LibShalom.h"

static double gtod_ref_time_sec = 0.0;

double dclock()
{
        double the_time, norm_sec;
        struct timeval tv;

        gettimeofday( &tv, NULL );

        if ( gtod_ref_time_sec == 0.0 )
                gtod_ref_time_sec = ( double ) tv.tv_sec;

        norm_sec = ( double ) tv.tv_sec - gtod_ref_time_sec;

        the_time = norm_sec + tv.tv_usec * 1.0e-6;

        return the_time;
}

void random_matrix( int m, int n, float *a)
{
  double drand48();
  int i,j;

  for ( i=0; i<m; i++ )
    for ( j=0; j<n; j++ )
		a[i*n+j]= 2.0 * (float)drand48( ) - 1.0 ;
}

void random_matrix1( int m, int n, float *a)
{
  double drand48();
  int i,j;

  for ( i=0; i<m; i++ )
    for ( j=0; j<n; j++ )
		a[i*n+j]= 1.0 ;
}

void transpose( int m, int n, float *a)
{
	float *temp_a = ( float * ) malloc( m* n * sizeof( float ) );

	int i, j;

	for( i =0 ;i< m; i++)
	{
		for( j= 0;j < n; j++)
		{
			temp_a[j * m + i] = a[i*n + j];
		}
	}

	for( i =0 ; i< n; i++)
	{
		for( j =0; j< m; j++)
		{
			a[ i * m+j] = temp_a[i*m+j];
		}
	}
}

int main()
{

	LibShalom_set_thread_nums(64);
	int i,j,k;
	int loop= 5 ;
	long M, N, K;
	float *C1;
	double start, cost;
	double gflops;
	int flag =0 ;
	char transa = 'N', transb = 'N';
	int pc;

	FILE *fp;

  	if( (fp=fopen("NSMM_small.txt","w")) == NULL )
  	{
    	puts("Fail to open file!");
    	exit(0);
  	}

    for( j =8; j < 9; j = j + 8)
    {

        M = 256 ;
        N = 4096;
        K = 5000;
        //M= N = K = j;
    	double ops = M *N *K * 1.0e-09 * 2;
    	fprintf(fp, "%d %d %d", M,N,K);

    	for(pc =0 ;pc < 1; pc++)
    	{
    		int flag =0 ;
    		int lda = K;
			int ldb = K;
			int ldc = N;

	    	float *A = ( float * ) malloc( K* M * sizeof( float ) );
	    	float *B = ( float * ) malloc( K* N * sizeof( float ) );
	    	float *C = ( float * ) malloc( M* N * sizeof( float ) );
	    	float *C1 = ( float * ) malloc( M* N * sizeof( float ) );
	    	float *SB= ( float * ) malloc( K* N * sizeof( float ) );

	    	random_matrix(M,K,A);
	    	random_matrix(K,N,B);

	    	transpose(K, N, B);

	    	for( i =0 ;i< 2; i++)
			cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, M, N, K, 1.0, A, lda, B, ldb, 0.0, C1, ldc);

			start = dclock();
		    for(i= 0; i<loop ;i++)
		    {

		     	cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, M, N, K, 1.0, A, lda, B, ldb, 0.0, C1, ldc);
		    }

		    for( i =0 ;i< 5; i++)
		    	LibShalom_sgemm_mp(NoTrans, Trans, C, A, B, M, N, K);

			start = dclock();
			for( i= 0; i< loop ;i++)
				LibShalom_sgemm_mp(NoTrans, Trans, C, A, B, M, N, K);
			cost =(dclock()-start)/loop;

			ops = M * N  * K * 1.0e-09 * 2;
			printf("N_SMM:  M= %d N=%d K=%d flops = %lf effic= %.3lf %\n", M, N, K,
					ops/cost, ops/cost/17.6 * 100/ 64);
			fprintf(fp, " %.3f", ops/cost);

		for( i= 0; i< M; i++)
		{
			for( j= 0 ;j<N;j++)
			{
				if((C[i*N+j]- C1[i*N+j] > 0.0001 || C[i*N+j]- C1[i*N+j] < -0.0001) )
				{
					printf("i = %d, j= %d\n",i ,j );
					printf("C= %lf , C1= %lf\n", C[i*N+j], C1[i*N+j]);
					flag =1;
					printf("error\n");
					break;
				}
			}
			if(flag ==1)
				break;
		}

		if(flag == 0)
			printf("结果正确\n");

		    free(A);
		    free(B);
		    free(C);
		    free(SB);
		}

		fprintf(fp, "\n");
	}

	fclose(fp);
    return 0;
}

/*

int main()
{

	LibShalom_thread_nums(32);
//	printf("%d\n", vec.size());
	int i,j,k;
	int loop= 100 ;
	long M, N, K;
	float *C1;
	double start, cost;
	double gflops;
	int flag =0 ;
	char transa = 'N', transb = 'N';
	int pc;

	FILE *fp;

	long Me[5]={5,13,13,23,26},
		Ne[5]={5,5,13,23,26},
		Ke[5]={5,13,13,23,13};

  	if( (fp=fopen("NSMM_small.txt","w")) == NULL )
  	{
    	puts("Fail to open file!");
    	exit(0);
  	}

    for( j = 0; j < 5; j++)
    {

        M = Me[j];
        N = Ne[j];
        K = Ke[j];
        //M= N = K = j;
    	double ops = M *N *K * 1.0e-09 * 2;
    	fprintf(fp, "%d %d %d", M,N,K);

    	for(pc =0 ;pc < 5; pc++)
    	{
	    	double *A = ( double * ) malloc( K* M * sizeof( double ) );
	    	double *B = ( double * ) malloc( K* N * sizeof( double ) );
	    	double *C = ( double * ) malloc( M* N * sizeof( double ) );
	    	double *SB= ( double * ) malloc( K* N * sizeof( double ) );

		    for( i =0 ;i< 5; i++)
		    	LibShalom_dgemm(NoTrans, NoTrans, C, A, B, M, N, K);

			start = dclock();
			for( i= 0; i< loop ;i++)
				LibShalom_dgemm(NoTrans, NoTrans, C, A, B, M, N, K);
			cost =(dclock()-start)/loop;

			ops = M * N  * K * 1.0e-09 * 2;
			printf("N_SMM:  M= %d N=%d K=%d flops = %lf effic= %.3lf %\n", M, N, K,
					ops/cost, ops/cost/8.8 * 100);
			fprintf(fp, " %.3f", ops/cost);

		    free(A);
		    free(B);
		    free(C);
		    free(SB);
		}

		fprintf(fp, "\n");
	}

	fclose(fp);
    return 0;
}
*/

max_stars_repo_name: AnonymousYWL/MYLIB

max_issues_repo_name: AnonymousYWL/LibShalom

max_forks_repo_name: AnonymousYWL/LibShalom