text:
/* multilarge_nlinear/fdf.c
 *
 * Copyright (C) 2015, 2016 Patrick Alken
 *
 * 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multilarge_nlinear.h>

gsl_multilarge_nlinear_workspace *
gsl_multilarge_nlinear_alloc (const gsl_multilarge_nlinear_type * T,
                              const gsl_multilarge_nlinear_parameters * params,
                              const size_t n, const size_t p)
{
  gsl_multilarge_nlinear_workspace * w;

  if (n < p)
    {
      GSL_ERROR_VAL ("insufficient data points, n < p", GSL_EINVAL, 0);
    }

  w = calloc (1, sizeof (gsl_multilarge_nlinear_workspace));
  if (w == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for workspace",
                     GSL_ENOMEM, 0);
    }

  w->n = n;
  w->p = p;
  w->type = T;
  w->fdf = NULL;
  w->niter = 0;
  w->params = *params;

  /* the cgst method uses its own built-in linear solver */
  if (w->params.trs == gsl_multilarge_nlinear_trs_cgst)
    {
      w->params.solver = gsl_multilarge_nlinear_solver_none;
    }

  w->x = gsl_vector_calloc (p);
  if (w->x == 0)
    {
      gsl_multilarge_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
    }

  w->f = gsl_vector_calloc (n);
  if (w->f == 0)
    {
      gsl_multilarge_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0);
    }

  w->dx = gsl_vector_calloc (p);
  if (w->dx == 0)
    {
      gsl_multilarge_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0);
    }

  w->g = gsl_vector_alloc (p);
  if (w->g == 0)
    {
      gsl_multilarge_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for g", GSL_ENOMEM, 0);
    }

  if (w->params.solver == gsl_multilarge_nlinear_solver_cholesky ||
      w->params.solver == gsl_multilarge_nlinear_solver_mcholesky)
    {
      w->JTJ = gsl_matrix_alloc (p, p);
      if (w->JTJ == 0)
        {
          gsl_multilarge_nlinear_free (w);
          GSL_ERROR_VAL ("failed to allocate space for JTJ", GSL_ENOMEM, 0);
        }
    }

  w->sqrt_wts_work = gsl_vector_calloc (n);
  if (w->sqrt_wts_work == 0)
    {
      gsl_multilarge_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for weights", GSL_ENOMEM, 0);
    }

  w->state = (T->alloc)(&(w->params), n, p);
  if (w->state == 0)
    {
      gsl_multilarge_nlinear_free (w);
      GSL_ERROR_VAL ("failed to allocate space for state", GSL_ENOMEM, 0);
    }

  return w;
}

void
gsl_multilarge_nlinear_free (gsl_multilarge_nlinear_workspace * w)
{
  RETURN_IF_NULL (w);

  if (w->state)
    (w->type->free) (w->state);

  if (w->dx)
    gsl_vector_free (w->dx);

  if (w->x)
    gsl_vector_free (w->x);

  if (w->f)
    gsl_vector_free (w->f);

  if (w->sqrt_wts_work)
    gsl_vector_free (w->sqrt_wts_work);

  if (w->g)
    gsl_vector_free (w->g);

  if (w->JTJ)
    gsl_matrix_free (w->JTJ);

  free (w);
}

gsl_multilarge_nlinear_parameters
gsl_multilarge_nlinear_default_parameters(void)
{
  gsl_multilarge_nlinear_parameters params;

  params.trs = gsl_multilarge_nlinear_trs_lm;
  params.scale = gsl_multilarge_nlinear_scale_more;
  params.solver = gsl_multilarge_nlinear_solver_cholesky;
  params.fdtype = GSL_MULTILARGE_NLINEAR_FWDIFF;
  params.factor_up = 3.0;
  params.factor_down = 2.0;
  params.avmax = 0.75;
  params.h_df = GSL_SQRT_DBL_EPSILON;
  params.h_fvv = 0.01;
  params.max_iter = 0;
  params.tol = 1.0e-6;

  return params;
}

int
gsl_multilarge_nlinear_init (const gsl_vector * x,
                             gsl_multilarge_nlinear_fdf * fdf,
                             gsl_multilarge_nlinear_workspace * w)
{
  return gsl_multilarge_nlinear_winit(x, NULL, fdf, w);
}

int
gsl_multilarge_nlinear_winit (const gsl_vector * x,
                              const gsl_vector * wts,
                              gsl_multilarge_nlinear_fdf * fdf,
                              gsl_multilarge_nlinear_workspace * w)
{
  const size_t n = w->f->size;

  if (n != fdf->n)
    {
      GSL_ERROR ("function size does not match workspace", GSL_EBADLEN);
    }
  else if (w->x->size != x->size)
    {
      GSL_ERROR ("vector length does not match workspace", GSL_EBADLEN);
    }
  else if (wts != NULL && n != wts->size)
    {
      GSL_ERROR ("weight vector length does not match workspace", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* initialize counters for function and Jacobian evaluations */
      fdf->nevalf = 0;
      fdf->nevaldfu = 0;
      fdf->nevaldf2 = 0;
      fdf->nevalfvv = 0;

      w->fdf = fdf;
      gsl_vector_memcpy(w->x, x);
      w->niter = 0;

      if (wts)
        {
          w->sqrt_wts = w->sqrt_wts_work;

          for (i = 0; i < n; ++i)
            {
              double wi = gsl_vector_get(wts, i);
              gsl_vector_set(w->sqrt_wts, i, sqrt(wi));
            }
        }
      else
        {
          w->sqrt_wts = NULL;
        }

      return (w->type->init) (w->state, w->sqrt_wts, w->fdf,
                              w->x, w->f, w->g, w->JTJ);
    }
}

int
gsl_multilarge_nlinear_iterate (gsl_multilarge_nlinear_workspace * w)
{
  int status =
    (w->type->iterate) (w->state, w->sqrt_wts, w->fdf,
                        w->x, w->f, w->g, w->JTJ, w->dx);

  w->niter++;

  return status;
}

double
gsl_multilarge_nlinear_avratio (const gsl_multilarge_nlinear_workspace * w)
{
  return (w->type->avratio) (w->state);
}

int
gsl_multilarge_nlinear_rcond (double * rcond, const gsl_multilarge_nlinear_workspace * w)
{
  int status = (w->type->rcond) (rcond, w->JTJ, w->state);
  return status;
}

int
gsl_multilarge_nlinear_covar (gsl_matrix * covar, gsl_multilarge_nlinear_workspace * w)
{
  if (covar->size1 != covar->size2)
    {
      GSL_ERROR ("covariance matrix must be square", GSL_ENOTSQR);
    }
  else if (covar->size1 != w->p)
    {
      GSL_ERROR ("covariance matrix does not match workspace", GSL_EBADLEN);
    }
  else
    {
      int status = (w->type->covar) (w->JTJ, covar, w->state);
      return status;
    }
}

/*
gsl_multilarge_nlinear_driver()
  Iterate the nonlinear least squares solver until completion

Inputs: maxiter  - maximum iterations to allow
        xtol     - tolerance in step x
        gtol     - tolerance in gradient
        ftol     - tolerance in ||f||
        callback - callback function to call each iteration
        callback_params - parameters to pass to callback function
        info     - (output) info flag on why iteration terminated
                   1 = stopped due to small step size ||dx|
                   2 = stopped due to small gradient
                   3 = stopped due to small change in f
                   GSL_ETOLX = ||dx|| has converged to within machine
                               precision (and xtol is too small)
                   GSL_ETOLG = ||g||_inf is smaller than machine
                               precision (gtol is too small)
                   GSL_ETOLF = change in ||f|| is smaller than machine
                               precision (ftol is too small)
        w        - workspace

Return:
GSL_SUCCESS if converged
GSL_MAXITER if maxiter exceeded without converging
GSL_ENOPROG if no accepted step found on first iteration
*/

int
gsl_multilarge_nlinear_driver (const size_t maxiter,
                               const double xtol,
                               const double gtol,
                               const double ftol,
                               void (*callback)(const size_t iter, void *params,
                                                const gsl_multilarge_nlinear_workspace *w),
                               void *callback_params,
                               int *info,
                               gsl_multilarge_nlinear_workspace * w)
{
  int status;
  size_t iter = 0;

  /* call user callback function prior to any iterations
   * with initial system state */
  if (callback)
    callback(iter, callback_params, w);

  do
    {
      status = gsl_multilarge_nlinear_iterate (w);

      /*
       * If the solver reports no progress on the first iteration,
       * then it didn't find a single step to reduce the
       * cost function and more iterations won't help so return.
       *
       * If we get a no progress flag on subsequent iterations,
       * it means we did find a good step in a previous iteration,
       * so continue iterating since the solver has now reset
       * mu to its initial value.
       */
      if (status == GSL_ENOPROG && iter == 0)
        {
          *info = status;
          return GSL_EMAXITER;
        }

      ++iter;

      if (callback)
        callback(iter, callback_params, w);

      /* test for convergence */
      status = gsl_multilarge_nlinear_test(xtol, gtol, ftol, info, w);
    }
  while (status == GSL_CONTINUE && iter < maxiter);

  /*
   * the following error codes mean that the solution has converged
   * to within machine precision, so record the error code in info
   * and return success
   */
  if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG)
    {
      *info = status;
      status = GSL_SUCCESS;
    }

  /* check if max iterations reached */
  if (iter >= maxiter && status != GSL_SUCCESS)
    status = GSL_EMAXITER;

  return status;
} /* gsl_multilarge_nlinear_driver() */

const char *
gsl_multilarge_nlinear_name (const gsl_multilarge_nlinear_workspace * w)
{
  return w->type->name;
}

gsl_vector *
gsl_multilarge_nlinear_position (const gsl_multilarge_nlinear_workspace * w)
{
  return w->x;
}

gsl_vector *
gsl_multilarge_nlinear_residual (const gsl_multilarge_nlinear_workspace * w)
{
  return w->f;
}

gsl_vector *
gsl_multilarge_nlinear_step (const gsl_multilarge_nlinear_workspace * w)
{
  return w->dx;
}

size_t
gsl_multilarge_nlinear_niter (const gsl_multilarge_nlinear_workspace * w)
{
  return w->niter;
}

const char *
gsl_multilarge_nlinear_trs_name (const gsl_multilarge_nlinear_workspace * w)
{
  return w->params.trs->name;
}

/*
gsl_multilarge_nlinear_eval_f()
  Compute residual vector y with user callback function, and apply
weighting transform if given:

y~ = sqrt(W) y

Inputs: fdf  - callback function
        x    - model parameters
        swts - weight matrix sqrt(W) = sqrt(diag(w1,w2,...,wn))
               set to NULL for unweighted fit
        y    - (output) (weighted) residual vector
               y_i = sqrt(w_i) f_i where f_i is unweighted residual
*/

int
gsl_multilarge_nlinear_eval_f(gsl_multilarge_nlinear_fdf *fdf,
                              const gsl_vector *x,
                              const gsl_vector *swts,
                              gsl_vector *y)
{
  int s = ((*((fdf)->f)) (x, fdf->params, y));

  ++(fdf->nevalf);

  /* y <- sqrt(W) y */
  if (swts)
    gsl_vector_mul(y, swts);

  return s;
}

/*
gsl_multilarge_nlinear_eval_df()
  Compute Jacobian matrix-vector product:

v = J * u

or

v = J^T u

Inputs: TransJ - use J or J^T
        x      - model parameters
        f      - residual vector f(x)
        u      - input vector u
        swts   - weight matrix W = diag(w1,w2,...,wn)
                 set to NULL for unweighted fit
        h      - finite difference step size
        fdtype - finite difference method
        fdf    - callback function
        v      - (output) vector v
        JTJ    - (output) matrix J^T J
        work   - workspace for finite difference, size n
*/

int
gsl_multilarge_nlinear_eval_df(const CBLAS_TRANSPOSE_t TransJ,
                               const gsl_vector *x,
                               const gsl_vector *f,
                               const gsl_vector *u,
                               const gsl_vector *swts,
                               const double h,
                               const gsl_multilarge_nlinear_fdtype fdtype,
                               gsl_multilarge_nlinear_fdf *fdf,
                               gsl_vector *v,
                               gsl_matrix *JTJ,
                               gsl_vector *work)
{
  const size_t n = fdf->n;
  const size_t p = fdf->p;

  if (u != NULL && ((TransJ == CblasNoTrans && u->size != p) ||
                    (TransJ == CblasTrans && u->size != n)))
    {
      GSL_ERROR("u vector has wrong size", GSL_EBADLEN);
    }
  else if (v != NULL && ((TransJ == CblasNoTrans && v->size != n) ||
                         (TransJ == CblasTrans && v->size != p)))
    {
      GSL_ERROR("v vector has wrong size", GSL_EBADLEN);
    }
  else if (JTJ != NULL && ((JTJ->size1 != p) || (JTJ->size2 != p)))
    {
      GSL_ERROR("JTJ matrix has wrong size", GSL_EBADLEN);
    }
  else
    {
      int status = GSL_SUCCESS;

      if (fdf->df)
        {
          /* call user-supplied function */
          status = ((*((fdf)->df)) (TransJ, x, u, fdf->params, v, JTJ));

          if (v)
            ++(fdf->nevaldfu);

          if (JTJ)
            ++(fdf->nevaldf2);
        }
      else
        {
#if 0
          /* use finite difference Jacobian approximation */
          status = gsl_multilarge_nlinear_df(h, fdtype, x, swts, fdf, f, df, work);
#endif
        }

      return status;
    }
}

/*
gsl_multilarge_nlinear_eval_fvv()
  Compute second direction derivative vector yvv with user
callback function, and apply weighting transform if given:

yvv~ = sqrt(W) yvv

Inputs: h    - step size for finite difference, if needed
        x    - model parameters, size p
        v    - unscaled geodesic velocity vector, size p
        f    - residual vector f(x), size n
        swts - weight matrix sqrt(W) = sqrt(diag(w1,w2,...,wn))
               set to NULL for unweighted fit
        fdf  - callback function
        yvv  - (output) (weighted) second directional derivative vector
               yvv_i = sqrt(w_i) fvv_i where f_i is unweighted
        work - workspace, size p
*/

int
gsl_multilarge_nlinear_eval_fvv(const double h,
                                const gsl_vector *x,
                                const gsl_vector *v,
                                const gsl_vector *f,
                                const gsl_vector *swts,
                                gsl_multilarge_nlinear_fdf *fdf,
                                gsl_vector *yvv,
                                gsl_vector *work)
{
  int status = GSL_SUCCESS;

  if (fdf->fvv != NULL)
    {
      /* call user-supplied function */
      status = ((*((fdf)->fvv)) (x, v, fdf->params, yvv));
      ++(fdf->nevalfvv);
    }
  else
    {
#if 0
      /* use finite difference approximation */
      status = gsl_multilarge_nlinear_fdfvv(h, x, v, f, J,
                                          swts, fdf, yvv, work);
#endif
    }

  /* yvv <- sqrt(W) yvv */
  if (swts)
    gsl_vector_mul(yvv, swts);

  return status;
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <inttypes.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <cblas.h>
#include <string.h> // memcpy
#include <math.h>
#include <clblast_c.h>

#include "cgraph.h"
#include "cg_operation.h"
#include "cg_types.h"
#include "cg_variables.h"
#include "cg_errors.h"
#include "cg_constants.h"
#include "cg_enums.h"
#include "cg_factory.h"
#include "cg_math.h"

//#pragma OPENCL EXTENSION cl_khr_fp64 : enable

#if cg_float == float
#define cblas_dcopy cblas_scopy
#define cblas_dscal cblas_sscal
#define cblas_dgemm cblas_sgemm
#define cblas_ddot cblas_sdot
#define cblas_dgemv cblas_sgemv
#define cl_double cl_float
#endif

void checkErr(CCLErr* err){
    if (err != NULL) {
        fprintf(stderr, "\n%s\n", err->message);
        exit(-1);
    }
}

#define CHECK_ERROR(err) \
    checkErr(err);
    //if (err != NULL) { fprintf(stderr, "\n%s\n", err->message); exit(-1); }

CCLContext * ctx = NULL;
CCLProgram* prg = NULL;
CCLQueue* queue = NULL;
CCLDevice* dev = NULL;

#if cg_float == float
#define CLBlastDcopy CLBlastScopy
#define CLBlastDscal CLBlastSscal
#define CLBlastDgemm CLBlastSgemm
#define CLBlastDdot CLBlastSdot
#define CLBlastDomatcopy CLBlastSomatcopy
#define CLBlastDsum CLBlastSsum
#endif

void copyDataToHost(CGResultNode* res){
    CCLEvent* evt = NULL;
    CCLErr * err = NULL;

    switch(res->type){
        case CGVT_DOUBLE:
            break;

        case CGVT_VECTOR:{
            CGVector* V = (CGVector*)res->value;
            V->data = malloc(V->len* sizeof(cgcl_float));
            if (V->loc == CG_DATALOC_DEVICE_MEM) {
                evt = ccl_buffer_enqueue_read(V->buf, queue, CL_TRUE, 0,
                                              V->len * sizeof(cgcl_float), V->data, NULL, NULL);
                if (!evt) exit(-1);
                ccl_buffer_destroy(V->buf);
            }
            else
                printf("FUCK\n");
            V->loc = CG_DATALOC_HOST_MEM;
            break;
        }

        case CGVT_MATRIX:{
            CGMatrix* M = (CGMatrix*)res->value;
            M->data = malloc(M->rows*M->cols* sizeof(cgcl_float));

            if (M->loc == CG_DATALOC_DEVICE_MEM) {
                evt = ccl_buffer_enqueue_read(M->buf, queue, CL_TRUE, 0,
                                              M->rows * M->cols * sizeof(cgcl_float), M->data, NULL, NULL);
                if (!evt) exit(-1);
                ccl_buffer_destroy(M->buf);
            }
            else
            printf("FUCK2\n");
            M->loc = CG_DATALOC_HOST_MEM;
            break;
        }
    }

    cl_bool status = ccl_queue_finish(queue, NULL);
    //if (!status) exit(-1);
    ////ccl_queue_destroy(queue);
}

CGDouble* makeDeviceDouble(){
    CCLErr * err = NULL;
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    cg_float* ptr = &Y->value;
    Y->buf = ccl_buffer_new(ctx, CL_MEM_READ_WRITE,
                            sizeof(cgcl_float), NULL, &err);
    CHECK_ERROR(err)

    Y->loc = CG_DATALOC_DEVICE_MEM;
    Y->value = 0.0;

    return Y;
}

CGVector* makeDeviceVector(uint64_t len){
    CCLErr * err = NULL;
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = len;
    Y->data = NULL; //calloc(len, sizeof(cg_float));

    Y->buf = ccl_buffer_new(ctx, CL_MEM_READ_WRITE ,
                            len * sizeof(cgcl_float), Y->data, &err);

    CHECK_ERROR(err)

    Y->loc = CG_DATALOC_DEVICE_MEM;

    return Y;
}

CGVector* makeDeviceVectorZeros(uint64_t len){
    CCLErr * err = NULL;
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = len;
    Y->data = calloc(len, sizeof(cg_float));

    Y->buf = ccl_buffer_new(ctx, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
                            len * sizeof(cgcl_float), Y->data, &err);

    CHECK_ERROR(err)

    Y->loc = CG_DATALOC_DEVICE_MEM;

    return Y;
}

CGMatrix* makeDeviceMatrix(uint64_t rows, uint64_t cols){
    CCLErr * err = NULL;
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = rows;
    Y->cols = cols;
    Y->data = NULL; //calloc(rows*cols, sizeof(cg_float));

    Y->buf = ccl_buffer_new(ctx, CL_MEM_READ_WRITE ,
                            rows*cols * sizeof(cgcl_float), Y->data, &err);
    CHECK_ERROR(err)

    Y->loc = CG_DATALOC_DEVICE_MEM;

    return Y;
}

CGMatrix* makeDeviceMatrixZeros(uint64_t rows, uint64_t cols){
    CCLErr * err = NULL;
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = rows;
    Y->cols = cols;
    Y->data = calloc(rows*cols, sizeof(cg_float));

    Y->buf = ccl_buffer_new(ctx, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
                            rows*cols * sizeof(cgcl_float), Y->data, &err);
    CHECK_ERROR(err)

    Y->loc = CG_DATALOC_DEVICE_MEM;

    return Y;
}

int selectContext() {
    /* Code. */
    ctx = ccl_context_new_from_menu(NULL);
    if (ctx == NULL)
        exit(-1);

    CCLErr * err = NULL;
    prg = ccl_program_new_from_source_file(ctx, "/Users/praisethemoon/projects/cgraph/source/libcgraph/kernels/vec.cl", NULL);
    cl_bool status = ccl_program_build(prg, NULL, &err);
    printf("%s\n", ccl_program_get_build_log(prg, NULL));
    CHECK_ERROR(err)

    dev = ccl_context_get_device(ctx, 0, &err);
    queue  = ccl_queue_new(ctx, dev, 0, &err);

    return 0;
}

/*
 * *********************
 * Multiplication
 * *********************
 */

/*
 * d.d
 */
CGResultNode* mulDD(CGDouble* M, CGDouble* V, CGraph* graph, CGNode* parentNode){
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = M->value * V->value;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * d.V == V.d
 */
CGResultNode* mulDV(CGDouble* A, CGVector* V, CGraph* graph, CGNode* parentNode){
    CGVector* Y = makeDeviceVector(V->len);
    CCLErr * err = NULL;

    //if(A->value != 0){
        uint64_t len = V->len;
        CCLBuffer * a = V->buf, * c = NULL;

        CHECK_ERROR(err)

        CHECK_ERROR(err)

        CCLEvent* evt = NULL;
        size_t gws = len;

        cl_mem V_mem = ccl_buffer_unwrap(V->buf);
        cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

        cl_command_queue q = ccl_queue_unwrap(queue);

        //cblas_dcopy(V->len, V->data, 1, y, 1);
        //cblas_dscal(V->len, a->value, y, 1);

        CLBlastDcopy(V->len, V_mem, 0.0, 1, Y_mem, 0.0, 1, &q, NULL);
        CLBlastDscal(V->len, A->value, Y_mem, 0.0, 1, &q, NULL);
        //ccl_queue_destroy(queue);
    //}

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * d.M == M.d
 */
CGResultNode* mulDM(CGDouble* A, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->cols*M->rows;
    cg_float value = A->value;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    //if(A->value != 0){
        CCLErr * err = NULL;

        CCLBuffer * a = M->buf;

        CHECK_ERROR(err)

        CHECK_ERROR(err)

        cl_mem M_mem = ccl_buffer_unwrap(M->buf);
        cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

        cl_command_queue q = ccl_queue_unwrap(queue);

        //cblas_dcopy(V->len, V->data, 1, y, 1);
        //cblas_dscal(V->len, a->value, y, 1);

        CLBlastDcopy(size, M_mem, 0.0, 1, Y_mem, 0.0, 1, &q, NULL);
        CLBlastDscal(size, A->value, Y_mem, 0.0, 1, &q, NULL);

        //ccl_queue_destroy(queue);
    //}

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M*v
 * broadcasting
 */
CGResultNode* mulMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->cols != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Add M(%"PRIu64", %"PRIu64") by V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    uint64_t i = 0;

    //res[i] = M->data[i] * V->data[(i)%V->len];

    uint64_t len = M->rows*M->cols;
    uint64_t vlen = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "mul_mv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, V->buf, Y->buf, ccl_arg_priv(vlen, cl_ulong), NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M*N
 * broadcasting
 */
CGResultNode* mulMM(CGMatrix* M1, CGMatrix* M2, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;
    if((M1->rows != M2->rows) && (M1->cols != M2->cols)){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot multiply M(%"PRIu64", %"PRIu64") by N(%"PRIu64", %"PRIu64")", M1->rows, M1->cols, M2->rows, M2->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M1->cols*M1->rows;

    CGMatrix* Y = makeDeviceMatrix(M1->rows, M1->cols);
    uint64_t len = M1->rows*M1->cols;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "mul_vv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M1->buf, M2->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M.v
 */
CGResultNode* dotMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;
    if(M->cols != 1){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot calculate M(%"PRIu64", %"PRIu64") DOT V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t len =  M->rows*V->len;

    CGMatrix* Y = makeDeviceMatrix(M->rows, V->len);

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    cl_command_queue q = ccl_queue_unwrap(queue);
    cl_mem M_mem = ccl_buffer_unwrap(M->buf);
    cl_mem V_mem = ccl_buffer_unwrap(V->buf);
    cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

    CLBlastDgemm(CLBlastLayoutRowMajor,
                 CLBlastTransposeNo, CLBlastTransposeNo, M->rows, V->len, M->cols,
                 1.0,
                 M_mem, 0, M->cols,
                 V_mem, 0, V->len,
                 0.0,
                 Y_mem, 0, Y->cols,
                 &q, NULL
                 );

    //evt = ccl_buffer_enqueue_read(Y->buf, queue, CL_TRUE, 0,
    //                              len * sizeof(cl_double), Y->data, NULL, NULL);

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * v.M
 */
CGResultNode* dotVM(CGVector* V, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;
    if(M->rows != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot calculate V(%"PRIu64") DOT M(%"PRIu64", %"PRIu64")", V->len, M->rows, M->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    CGVector* Y = makeDeviceVector(M->cols);

    uint64_t m = 1;
    uint64_t n = M->cols;
    uint64_t k = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    cl_command_queue q = ccl_queue_unwrap(queue);
    cl_mem M_mem = ccl_buffer_unwrap(M->buf);
    cl_mem V_mem = ccl_buffer_unwrap(V->buf);
    cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

    //cblas_dgemm(CblasRowMajor,
    //            CblasNoTrans, CblasNoTrans, m, n, k, 1, V->data, k, M->data, n, 0, Y->data, n);

    CLBlastDgemm(CLBlastLayoutRowMajor,
                 CLBlastTransposeNo, CLBlastTransposeNo, m, n, k,
                 1.0,
                 V_mem, 0, k,
                 M_mem, 0, n,
                 0.0,
                 Y_mem, 0, n,
                 &q, NULL
    );

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M.N
 */
CGResultNode* dotMM(CGMatrix* M, CGMatrix* N, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;
    if(M->cols != N->rows){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot DOT M(%"PRIu64", %"PRIu64") by N(%"PRIu64", %"PRIu64")", M->rows, M->cols, N->rows, N->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    CGMatrix* Z = makeDeviceMatrix(M->rows, N->cols);

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    cl_command_queue q = ccl_queue_unwrap(queue);
    cl_mem M_mem = ccl_buffer_unwrap(M->buf);
    cl_mem N_mem = ccl_buffer_unwrap(N->buf);
    cl_mem Z_mem = ccl_buffer_unwrap(Z->buf);

    /*
    cblas_dgemm(CblasRowMajor,
                CblasNoTrans, CblasNoTrans, M->rows, N->cols, M->cols,
                1.0, M->data, M->cols, N->data, N->cols, 0, z, Z->cols);
    */

    CLBlastDgemm(CLBlastLayoutRowMajor,
                 CLBlastTransposeNo, CLBlastTransposeNo, M->rows, N->cols, M->cols,
                 1.0,
                 M_mem, 0, M->cols,
                 N_mem, 0, N->cols,
                 0.0,
                 Z_mem, 0, Z->cols,
                 &q, NULL
    );

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Z;

    return result;
}

/*
 * V.V
 */
CGResultNode* crossVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot compute Cross Product of V(%"PRIu64") by V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = V1->len;

    CGVector* Y = makeDeviceVector(V1->len);
    uint64_t len = V1->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "mul_vv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V1->buf, V2->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Dot Product
 * *********************
 */

/*
 * V.V
 */
CGResultNode* dotVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot compute Dot Product of V(%"PRIu64") by V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = V1->len;
    cg_float res = 0;

    CGDouble* Y = makeDeviceDouble();

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    cl_command_queue q = ccl_queue_unwrap(queue);
    cl_mem M_mem = ccl_buffer_unwrap(V1->buf);
    cl_mem N_mem = ccl_buffer_unwrap(V2->buf);
    cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

    //Y->value = cblas_ddot(V1->len, V1->data, 1, V2->data, 1);

    CLBlastDdot(size, Y_mem, 0,
            M_mem, 0.0, 1,
            N_mem, 0.0, 1,
            &q, NULL);

    evt = ccl_buffer_enqueue_read(Y->buf, queue, CL_TRUE, 0,
                                  sizeof(cl_double), &Y->value, NULL, NULL);

    //ccl_queue_destroy(queue);
    ccl_buffer_destroy(Y->buf);

    Y->loc = CG_DATALOC_HOST_MEM;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Division
 * *********************
 */

/*
 *  d/d
 */
CGResultNode* divDD(CGDouble* D1, CGDouble* D2, CGraph* graph, CGNode* parentNode){
    if(D2->value == 0.0){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Diving by zero D/(D=0)");
        return returnResultError(graph, CGET_DIVIDE_BY_ZERO, parentNode, msg);
    }

    cg_float res = D1->value / D2->value;

    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = res;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 *  V/d
 */
CGResultNode* divVD(CGVector* V, CGDouble* D, CGraph* graph, CGNode* parentNode){
    if(D->value == 0.0){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Diving by zero V/(D=0)");
        return returnResultError(graph, CGET_DIVIDE_BY_ZERO, parentNode, msg);
    }

    CGVector* Y = makeDeviceVector(V->len);
    CCLErr * err = NULL;

    //if(A->value != 0){
    uint64_t len = V->len;
    CCLBuffer * a = V->buf, * c = NULL;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    cl_mem V_mem = ccl_buffer_unwrap(V->buf);
    cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

    cl_command_queue q = ccl_queue_unwrap(queue);

    //cblas_dcopy(V->len, V->data, 1, y, 1);
    //cblas_dscal(V->len, a->value, y, 1);

    CLBlastDcopy(V->len, V_mem, 0.0, 1, Y_mem, 0.0, 1, &q, NULL);
    CLBlastDscal(V->len, 1.0/D->value, Y_mem, 0.0, 1, &q, NULL);
    //ccl_queue_destroy(queue);
    //}

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;

}

/*
 *  V/V
 */
CGResultNode* divVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot calculate  V(%"PRIu64") DIV V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }
    CCLErr * err = NULL;

    uint64_t size = V1->len;

    CGVector* Y = makeDeviceVector(V1->len);
    uint64_t len = V1->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "div_vv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V1->buf, V2->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 *  d/V element-wise
 */
CGResultNode* divDV(CGDouble* D, CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    uint64_t size = V->len;

    CGVector* Y = makeDeviceVector(V->len);
    uint64_t len = V->len;

    // TODO: check if D.value == 0.0

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "div_dv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), ccl_arg_priv(D->value, cgcl_float), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 *  M/d
 */
CGResultNode* divMD(CGMatrix* M, CGDouble* D, CGraph* graph, CGNode* parentNode){
    if(D->value == 0.0){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Diving by zero M/(D=0)");
        return returnResultError(graph, CGET_DIVIDE_BY_ZERO, parentNode, msg);
    }

    uint64_t size = M->cols*M->rows;
    cg_float value = D->value;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    //if(A->value != 0){
    CCLErr * err = NULL;

    CCLBuffer * a = M->buf;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    cl_mem M_mem = ccl_buffer_unwrap(M->buf);
    cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

    cl_command_queue q = ccl_queue_unwrap(queue);

    CLBlastDcopy(size, M_mem, 0.0, 1, Y_mem, 0.0, 1, &q, NULL);
    CLBlastDscal(size, 1.0/D->value, Y_mem, 0.0, 1, &q, NULL);

    //ccl_queue_destroy(queue);
    //}

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M/v
 */
CGResultNode* divMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->rows != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot calculate  M(%"PRIu64", %"PRIu64") DIV V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    uint64_t len = M->rows*M->cols;
    uint64_t vlen = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "div_mv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, V->buf, Y->buf, ccl_arg_priv(vlen, cl_ulong), NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 *  d/M element-wise
 */
CGResultNode* divDM(CGDouble* D, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t len = M->rows*M->cols;
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    // TODO: check if D.value == 0.0

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "div_dv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), ccl_arg_priv(D->value, cgcl_float), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Addition
 * *********************
 */

/*
 * d+d
 */
CGResultNode* addDD(CGDouble* D1, CGDouble* D2, CGraph* graph, CGNode* parentNode){
    cg_float res = D1->value + D2->value;

    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = res;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * V+d
 */
CGResultNode* addVD(CGVector* V, CGDouble* D, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    uint64_t size = V->len;

    CGVector* Y = makeDeviceVector(V->len);
    uint64_t len = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "add_dv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), ccl_arg_priv(D->value, cgcl_float), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M+d
 */
CGResultNode* addMD(CGMatrix* M, CGDouble* D, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    uint64_t size = M->rows*M->cols;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = size;

    evt = ccl_program_enqueue_kernel(prg, "add_dv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(size, cl_uint), ccl_arg_priv(D->value, cgcl_float), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;;
}

/*
 * V+V
 */
CGResultNode* addVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Add V(%"PRIu64") by V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = V1->len;
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(size);

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = size;

    evt = ccl_program_enqueue_kernel(prg, "add_vv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(size, cl_uint), V1->buf, V2->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M+V
 */
CGResultNode* addMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->cols != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Add M(%"PRIu64", %"PRIu64") by V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    uint64_t len = M->rows*M->cols;
    uint64_t vlen = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "add_mv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, V->buf, Y->buf, ccl_arg_priv(vlen, cl_ulong), NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M+M
 */
CGResultNode* addMM(CGMatrix* M1, CGMatrix* M2, CGraph* graph, CGNode* parentNode){
    if((M1->rows != M2->rows) && (M1->cols != M2->cols)){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot add M(%"PRIu64", %"PRIu64") by N(%"PRIu64", %"PRIu64")", M1->rows, M1->cols, M2->rows, M2->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    CCLErr * err = NULL;
    uint64_t size = M1->cols*M1->rows;

    CGMatrix* Y = makeDeviceMatrix(M1->rows, M1->cols);
    uint64_t len = M1->rows*M1->cols;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "add_vv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M1->buf, M2->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Substraction
 * *********************
 */

/*
 * d-d
 */
CGResultNode* subDD(CGDouble* D1, CGDouble* D2, CGraph* graph, CGNode* parentNode){
    cg_float res = D1->value - D2->value;

    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = res;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * V-d
 */
CGResultNode* subVD(CGVector* V, CGDouble* D, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    uint64_t size = V->len;

    CGVector* Y = makeDeviceVector(V->len);
    uint64_t len = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "sub_vd", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), ccl_arg_priv(D->value, cgcl_float), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * d-V
 */
CGResultNode* subDV(CGDouble* D, CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    uint64_t size = V->len;

    CGVector* Y = makeDeviceVector(V->len);
    uint64_t len = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "sub_dv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), ccl_arg_priv(D->value, cgcl_float), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M-d
 */
CGResultNode* subMD(CGMatrix* M, CGDouble* D, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    uint64_t size = M->rows*M->cols;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = size;

    evt = ccl_program_enqueue_kernel(prg, "sub_vd", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(size, cl_uint), ccl_arg_priv(D->value, cgcl_float), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * D-M
 */
CGResultNode* subDM(CGDouble* D, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    uint64_t size = M->rows*M->cols;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = size;

    evt = ccl_program_enqueue_kernel(prg, "sub_dv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(size, cl_uint), ccl_arg_priv(D->value, cgcl_float), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * V-V
 */
CGResultNode* subVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Substract V(%"PRIu64") by V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = V1->len;
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(size);

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = size;

    evt = ccl_program_enqueue_kernel(prg, "sub_vv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(size, cl_uint), V1->buf, V2->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M-M
 */
CGResultNode* subMM(CGMatrix* M1, CGMatrix* M2, CGraph* graph, CGNode* parentNode){
    if((M1->rows != M2->rows) && (M1->cols != M2->cols)){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Substract M(%"PRIu64", %"PRIu64") by N(%"PRIu64", %"PRIu64")", M1->rows, M1->cols, M2->rows, M2->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    CCLErr * err = NULL;
    uint64_t size = M1->cols*M1->rows;

    CGMatrix* Y = makeDeviceMatrix(M1->rows, M1->cols);
    uint64_t len = M1->rows*M1->cols;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "sub_vv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M1->buf, M2->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M-V
 */
CGResultNode* subMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->rows != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Substract M(%"PRIu64", %"PRIu64") by V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    uint64_t len = M->rows*M->cols;
    uint64_t vlen = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "sub_mv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, V->buf, Y->buf, ccl_arg_priv(vlen, cl_ulong), NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * V-M
 */
CGResultNode* subVM(CGVector* V, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    if(M->rows != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Substract M(%"PRIu64", %"PRIu64") by V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    uint64_t len = M->rows*M->cols;
    uint64_t vlen = V->len;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "sub_mv", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, V->buf, Y->buf, ccl_arg_priv(vlen, cl_ulong), NULL);
    CHECK_ERROR(err)

    //ccl_queue_destroy(queue);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Power
 * *********************
 */

/*
 * d^d
 */
CGResultNode* powDD(CGDouble* D1, CGDouble* D2, CGraph* graph, CGNode* parentNode){
    cg_float res = pow(D1->value, D2->value);

    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = res;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * V^d
 */
CGResultNode* powVD(CGVector* V, CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float value = D->value;

    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "pow_vd", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), ccl_arg_priv(value, cgcl_float), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M^d
 */
CGResultNode* powMD(CGMatrix* M, CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float value = D->value;

    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);

    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "pow_vd", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), ccl_arg_priv(value, cgcl_float), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Transpose Multiplication
 * *********************
 */
/*
 * M^T.v
 * TODO
 */
CGResultNode* mulMtV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = V->len;
    Y->data = y;

    cblas_dgemv(CblasRowMajor,
                CblasTrans, M->rows, M->cols, 1.0, M->data,
                V->len, V->data, 1, 0.0, Y->data, 1);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * exp(d)
 */
CGResultNode* expD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = exp(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * exp(v)
 */
CGResultNode* expV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "exp_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * exp(M)
 */
CGResultNode* expM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);
    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "exp_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * log(d)
 */
CGResultNode* logD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = log(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * log(V)
 */
CGResultNode* logV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "log_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * exp(M)
 */
CGResultNode* logM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);
    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "log_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Sin
 * *********************
 */

/*
 * sin(d)
 */
CGResultNode* sinD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = sin(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * sin(V)
 */
CGResultNode* sinV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "sin_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * sin(M)
 */
CGResultNode* sinM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);
    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "sin_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * cos
 * *********************
 */

/*
 * cos(d)
 */
CGResultNode* cosD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = cos(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * cos(V)
 */
CGResultNode* cosV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "cos_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * sin(M)
 */
CGResultNode* cosM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);
    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "cos_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * tan
 * *********************
 */

/*
 * tan(d)
 */
CGResultNode* tanD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = tan(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * tan(V)
 */
CGResultNode* tanV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "tan_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * tan(M)
 */
CGResultNode* tanM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);
    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "tan_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * tanh
 * *********************
 */

/*
 * tanh(d)
 */
CGResultNode* tanhD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = tanh(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * tanh(V)
 */
CGResultNode* tanhV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "tanh_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * tanh(M)
 */
CGResultNode* tanhM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);
    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "tanh_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Invert
 * *********************
 */

// TODO

/*
 * *********************
 * Transpose
 * *********************
 */

/*
 * D^t
 */
CGResultNode* transposeD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    CGDouble* R = calloc(1, sizeof(CGDouble));
    R->value = D->value;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = R;

    return result;
}

/*
 * V^t
 */
CGResultNode* transposeV(CGVector* V, CGraph* graph, CGNode* parentNode){
    uint64_t size = V->len;

    CCLErr * err = NULL;
    CCLEvent* evt = NULL;
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = V->len;
    Y->cols = 1;
    Y->data = NULL;

    Y->buf = ccl_buffer_new(ctx, CL_MEM_READ_WRITE ,
                            V->len * sizeof(cgcl_float), NULL, &err);

    CHECK_ERROR(err)

    CCLEventWaitList ewl = NULL;
    evt = ccl_buffer_enqueue_copy(V->buf, Y->buf, queue, 0, 0, V->len * sizeof(cgcl_float), NULL, &err);

    CHECK_ERROR(err)
    ccl_event_wait(ccl_ewl(&ewl, evt, NULL), &err);
    CHECK_ERROR(err)

    Y->loc = CG_DATALOC_DEVICE_MEM;
    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M^t
 */
CGResultNode* transposeM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->cols*M->rows;
    CGMatrix* Y = makeDeviceMatrix(M->cols, M->rows);

    //if(A->value != 0){
    CCLErr * err = NULL;

    CCLBuffer * a = M->buf;

    CHECK_ERROR(err)

    CHECK_ERROR(err)

    cl_mem M_mem = ccl_buffer_unwrap(M->buf);
    cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

    cl_command_queue q = ccl_queue_unwrap(queue);

    CLBlastDomatcopy(CLBlastLayoutRowMajor, CLBlastTransposeYes, M->rows, M->cols, 1.0f, M_mem, 0, M->cols, Y_mem, 0, Y->cols, &q, NULL);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * SUM
 * *********************
 */

/*
 * sum(D)
 */
CGResultNode* sumD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    return makeDoubleResultNode(D->value);
}

/*
 * sum(V)
 */
CGResultNode* sumV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CGDouble* Y = makeDeviceDouble();
    CCLEvent* evt = NULL;
    cl_command_queue q = ccl_queue_unwrap(queue);
    cl_mem V_mem = ccl_buffer_unwrap(V->buf);
    cl_mem Y_mem = ccl_buffer_unwrap(Y->buf);

    uint64_t  size = V->len;

    CLBlastDsum(size, Y_mem, 0, V_mem, 0, 1, &q, NULL);

    evt = ccl_buffer_enqueue_read(Y->buf, queue, CL_TRUE, 0,
                                  sizeof(cl_double), &Y->value, NULL, NULL);

    //ccl_queue_destroy(queue);
    ccl_buffer_destroy(Y->buf);

    Y->loc = CG_DATALOC_HOST_MEM;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * TODO:
 * sum(M, axis=axis)
 */

CGResultNode* sumM(CGMatrix* M, CGraph* graph, CGNode* parentNode, uint8_t axis){
    CCLErr * err = NULL;
    CGVector* Y = NULL;
    CCLEvent* evt = NULL;
    uint64_t  size = M->rows*M->cols;
    if(axis == 0){
        uint64_t len = M->cols;
        Y = makeDeviceVectorZeros(len);

        size_t gws = size;

        evt = ccl_program_enqueue_kernel(prg, "sum_matrix_rows", queue, 1, NULL, &gws,
                                         NULL, NULL, &err, ccl_arg_priv(size, cl_uint), ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);

        CHECK_ERROR(err)
    }
    else {
        uint64_t len = M->rows;
        Y = makeDeviceVectorZeros(len);

        size_t gws = size;

        evt = ccl_program_enqueue_kernel(prg, "sum_matrix_cols", queue, 1, NULL, &gws,
                                         NULL, NULL, &err, ccl_arg_priv(size, cl_uint), ccl_arg_priv(M->cols, cl_uint), M->buf, Y->buf, NULL);

        CHECK_ERROR(err)
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

CGResultNode* max(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = ((CGDouble*)res->value)->value;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 1;
            cg_float m = v->data[0];

            for (;i < v->len; i++){
                if(v->data[i] > m)
                    m = v->data[i];
            }

            Y->value = m;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->cols;i++){
                    y[i] = M->data[i];
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(y[i%len] < M->data[i])
                        y[i%len] = M->data[i];
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows;i++){
                    y[i] = M->data[i*M->cols];
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(y[i/M->cols] < M->data[i])
                        y[i/M->cols] = M->data[i];
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

CGResultNode* min(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = ((CGDouble*)res->value)->value;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 1;
            cg_float m = v->data[0];

            for (;i < v->len; i++){
                if(v->data[i] > m)
                    m = v->data[i];
            }

            Y->value = m;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->cols;i++){
                    y[i] = M->data[i];
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(y[i%len] > M->data[i])
                        y[i%len] = M->data[i];
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows;i++){
                    y[i] = M->data[i*M->cols];
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(y[i/M->cols] > M->data[i])
                        y[i/M->cols] = M->data[i];
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

CGResultNode* argmax(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = ((CGDouble*)res->value)->value;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 1;
            uint64_t maxIdx = 0;

            for (;i < v->len; i++){
                if(v->data[i] > v->data[maxIdx])
                    maxIdx = i;
            }

            Y->value = maxIdx;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->cols;i++){
                    y[i] = i;
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(M->data[(uint64_t)y[i%len]] < M->data[i])
                        y[i%len] = i;
                }

                // transform vector representation into matrix
                for(i=0;i<M->cols;i++){
                    y[i] = ((int)y[i])/M->cols;
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows;i++){
                    y[i] = i*M->cols;
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(M->data[(uint64_t)y[i/M->cols]] < M->data[i])
                        y[i/M->cols] = i;
                }

                // transform vector representation into matrix
                for(i=0;i<M->rows;i++){
                    y[i] = ((int)y[i])%M->cols;
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

CGResultNode* argmin(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = 0;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 1;
            uint64_t minIdx = 0;

            for (;i < v->len; i++){
                if(v->data[i] > v->data[minIdx])
                    minIdx = (cg_float)i;
            }

            Y->value = minIdx;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->cols;i++){
                    y[i] = i;
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(M->data[(uint64_t)y[i%len]] > M->data[i])
                        y[i%len] = (cg_float)i;
                }

                // transform vector representation into matrix
                for(i=0;i<M->cols;i++){
                    y[i] = ((int)y[i])/M->cols;
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows;i++){
                    y[i] = i*M->cols;
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(M->data[(uint64_t)y[i/M->cols]] > M->data[i])
                        y[i/M->cols] = (cg_float)i;
                }

                // transform vector representation into matrix
                for(i=0;i<M->rows;i++){
                    y[i] = ((int)y[i])%M->cols;
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

CGResultNode* mean(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = ((CGDouble*)res->value)->value;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 0;
            cg_float m = v->data[0];

            for (;i < v->len; i++){
                m += v->data[i];
            }

            m /= v->len;

            Y->value = m;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows*M->cols;i++){
                    y[i%len] += M->data[i];
                }

                for(i=0;i<M->cols;i++){
                    y[i] /= M->rows;
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows*M->cols;i++){
                    y[i/M->cols] += M->data[i];
                }

                for(i=0;i<M->rows;i++){
                    y[i] /= M->cols;
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

/*
 * relu(d)
 */
CGResultNode* reluD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float val = D->value;
    return makeDoubleResultNode(__cg_relu(val));
}

/*
 * relu(V)
 */
CGResultNode* reluV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "relu_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * relu(M)
 */
CGResultNode* reluM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);
    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "relu_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * softplus(d)
 */
CGResultNode* softplusD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float val = D->value;
    return makeDoubleResultNode(__cg_softplus(val));
}

/*
 * softplus(V)
 */
CGResultNode* softplusV(CGVector* V, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGVector* Y = makeDeviceVector(V->len);

    uint64_t len = V->len;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "softplus_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), V->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * softplus(M)
 */
CGResultNode* softplusM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    CCLErr * err = NULL;

    CGMatrix* Y = makeDeviceMatrix(M->rows, M->cols);
    uint64_t len = M->rows*M->cols;

    CCLEvent* evt = NULL;
    size_t gws = len;

    evt = ccl_program_enqueue_kernel(prg, "softplus_v", queue, 1, NULL, &gws,
                                     NULL, NULL, &err, ccl_arg_priv(len, cl_uint), M->buf, Y->buf, NULL);
    CHECK_ERROR(err)

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

#if cg_float == float
#undef cblas_dcopy
#undef cblas_dscal
#undef cblas_dgemm
#undef cblas_ddot
#undef cblas_dgemv
#undef cl_double
#endif

#if cg_float == float
#undef CLBlastDcopy
#undef CLBlastDscal
#undef CLBlastDgemm
#undef CLBlastDdot
#undef CLBlastDsum
#endif

ext: c

max_stars_repo_name: Enehcruon/cgraph

max_issues_repo_name: Enehcruon/cgraph

max_forks_repo_name: praisethemoon/ccgraph


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch

 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, or (at your option)
 any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 02110-1301, USA.

*/

#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_sf_erf.h>
#include <math.h>
#include <stdlib.h>

/* ------------------------------------------------------ */
gsl_complex gsl_complex_step_real (gsl_complex a)
{
  gsl_complex z;

  if (GSL_REAL(a) < 0)
    {
      GSL_SET_COMPLEX (&z, 0, 0);
    }
  else
    {
      GSL_SET_COMPLEX (&z, 1, 0);
    }

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_min_real (gsl_complex a, gsl_complex b)
{
  gsl_complex z;
  double min;

  /* just consider real parts */
  min = GSL_REAL(a) < GSL_REAL(b) ? GSL_REAL(a) : GSL_REAL(b);
  GSL_SET_COMPLEX (&z, min, 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_max_real (gsl_complex a, gsl_complex b)
{
  gsl_complex z;
  double max;

  /* just consider real parts */
  max = GSL_REAL(a) > GSL_REAL(b) ? GSL_REAL(a) : GSL_REAL(b);
  GSL_SET_COMPLEX (&z, max, 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_carg (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, gsl_complex_arg(a), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_cabs (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, gsl_complex_abs(a), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_cabs2 (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, gsl_complex_abs2(a), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_clogabs (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, gsl_complex_logabs(a), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_erf (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, gsl_sf_erf(GSL_REAL(a)), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_arctan2 (gsl_complex a, gsl_complex b)
{
  gsl_complex z, p;

  if(GSL_REAL(b) != 0.0)
    {
      z = gsl_complex_arctan(gsl_complex_div(a, b));
      if(GSL_REAL(b) < 0.0){
	GSL_SET_COMPLEX (&p, M_PI, 0);
	if(GSL_REAL(a) >= 0.0)
	  z = gsl_complex_add(z, p);
	else
	  z = gsl_complex_sub(z, p);
      }
    }
  else
    {
      if(GSL_REAL(a) >= 0.0)
	{
	  GSL_SET_COMPLEX (&z, M_PI/2.0, 0.0);
	}
      else
	{
	  GSL_SET_COMPLEX (&z, -M_PI/2.0, 0.0);
	}
    }

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_realpart (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, GSL_REAL(a), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_imagpart (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, GSL_IMAG(a), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_round (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, trunc(GSL_REAL(a)), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_ceiling (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, ceil(GSL_REAL(a)), 0);

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_floor (gsl_complex a)
{
  gsl_complex z;

  GSL_SET_COMPLEX (&z, floor(GSL_REAL(a)), 0);

  return z;
}

/* ------------------------------------------------------ */
void gsl_complex_rand_seed(long a)
{
  srandom(a);
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_rand (gsl_complex a)
{

  double r = random()/((double) RAND_MAX);
  return gsl_complex_rect(r, 0.0);
}

max_stars_repo_name: shunsuke-sato/octopus

max_issues_repo_name: shunsuke-sato/octopus

max_forks_repo_name: shunsuke-sato/octopus


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef xsd_284c8467_1f5e_4ae4_b645_077382bf6ab4_h
#define xsd_284c8467_1f5e_4ae4_b645_077382bf6ab4_h

#include <gslib/config.h>
#include <gslib/xml.h>
#include <gslib/std.h>
#include <gslib/xsdtypes.h>

__gslib_begin__

//#define xsd_def_

class xsd_schema;
struct xsd_extern_schema;
class xsd_node;
class xsd_element;
class xsd_attribute;
class xsd_group;
class xsd_attribute_group;
class xsd_simple_type;
class xsd_complex_type;

typedef xml_wrapper xsd_wrapper;
typedef vector<xsd_wrapper*> xsd_wrappers;
typedef xmltree::iterator xsd_iterator;
typedef xmltree::const_iterator xsd_const_iterator;

template<class _ty>
struct xsd_schema_ns_hash:
    public std::unary_function<_ty, size_t>
{
    size_t operator()(_ty t) const
    {
        assert(t);
        return string_hash(t->get_namespace());
    }
};

template<class _ty>
struct xsd_schema_ns_equalto:
    public std::binary_function<_ty, _ty, bool>
{
    bool operator()(_ty t1, _ty t2) const
    {
        assert(t1 && t2);
        return t1->get_namespace() == t2->get_namespace();
    }
};

template<class _ty>
struct xsd_schema_loc_hash:
    public std::unary_function<_ty, size_t>
{
    size_t operator()(_ty t) const
    {
        assert(t);
        return string_hash(t->get_location());
    }
};

template<class _ty>
struct xsd_schema_loc_equalto:
    public std::binary_function<_ty, _ty, bool>
{
    bool operator()(_ty t1, _ty t2) const
    {
        assert(t1 && t2);
        return t1->get_location() == t2->get_location();
    }
};

template<class _ty>
struct xsd_name_hash:
    public std::unary_function<_ty, size_t>
{
    size_t operator()(_ty t) const
    {
        assert(t);
        return string_hash(t->get_name());
    }
};

template<class _ty>
struct xsd_name_equalto:
    public std::binary_function<_ty, _ty, bool>
{
    bool operator()(_ty t1, _ty t2) const
    {
        assert(t1 && t2);
        return t1->get_name() == t2->get_name();
    }
};

typedef vector<xsd_schema*> xsd_schemas;
typedef unordered_multiset<xsd_schema*, xsd_schema_ns_hash<xsd_schema*>, xsd_schema_ns_equalto<xsd_schema*> > xsd_schema_ns_map;
typedef unordered_set<xsd_schema*, xsd_schema_loc_hash<xsd_schema*>, xsd_schema_loc_equalto<xsd_schema*> > xsd_schema_loc_map;
typedef list<xsd_extern_schema> xsd_extern_schemas;
typedef unordered_set<xsd_node*, xsd_name_hash<xsd_node*>, xsd_name_equalto<xsd_node*> > xsd_node_map;
typedef unordered_multiset<xsd_node*, xsd_name_hash<xsd_node*>, xsd_name_equalto<xsd_node*> > xsd_node_multimap;

struct xsd_global
{
    xsd_schemas             packs;
    xsd_schema_ns_map       ns_map;
    xsd_schema_loc_map      loc_map;

public:
    ~xsd_global();

};

struct xsd_cg_config {};

extern xsd_schema* xsd_parse_schema(xsd_global& global, xmltree& dom, const gchar* location);
extern xsd_schema* xsd_parse_schema(xsd_global& global, const gchar* location);
extern void xsd_make_file_name(string& name, const string& str);
extern bool xsd_prepare_name_mangling(xsd_global& global);
extern bool xsd_prepare_translation(xsd_schema* schema, xsd_global& global);
extern bool xsd_translate_declarations(string& output, xsd_schema* schema, xsd_global& global);
extern bool xsd_translate_implementations(string& output, xsd_schema* schema, xsd_global& global);

enum xsd_extern_schema_type
{
    xest_include,
    xest_import,
};

enum xsd_node_type
{
    xnt_key,
    xnt_element,
    xnt_attribute,
    xnt_attribute_group,
    xnt_group,
    xnt_simple_type,
    xnt_complex_type,
};

struct xsd_extern_schema
{
    xsd_extern_schema_type  type;
    string                  ns;
    xsd_schema*             schema;

public:
    xsd_schema* operator->() { return schema; }
    const xsd_schema* operator->() const { return schema; }
};

class xsd_schema
{
protected:
    xsd_wrapper*            _source;
    string                  _namespace;
    string                  _location;
    xsd_extern_schemas      _external;      // include & import
    string                  _xsd_namespace;
    xsd_node_map            _element_map;
    xsd_node_map            _attribute_map;
    xsd_node_map            _attribute_group_map;
    xsd_node_map            _group_map;
    xsd_node_map            _simple_type_map;
    xsd_node_map            _complex_type_map;

public:
    xsd_schema();
    ~xsd_schema();
    const string& get_namespace() const { return _namespace; }
    const string& get_location() const { return _location; }
    const string& get_xsd_namespace() const { return _xsd_namespace; }
    void set_namespace(const string& str) { _namespace = str; }
    void set_location(const string& str) { _location = str; }
    void set_xsd_namespace(const string& str) { _xsd_namespace = str; }
    void set_source(xsd_wrapper* s) { _source = s; }
    xsd_extern_schemas& get_extern_schemas() { return _external; }
    xsd_node_map& get_element_map() { return _element_map; }
    xsd_node_map& get_attribute_map() { return _attribute_map; }
    xsd_node_map& get_attribute_group_map() { return _attribute_group_map; }
    xsd_node_map& get_group_map() { return _group_map; }
    xsd_node_map& get_simple_type_map() { return _simple_type_map; }
    xsd_node_map& get_complex_type_map() { return _complex_type_map; }

protected:
    static xsd_node* find_node(const xsd_node_map& m, const string& name);
    template<class _node>
    static _node* find_and_verify(const xsd_node_map& m, const string& n, xsd_node_type t)
    {
        xsd_node* r = find_node(m, n);
        if(!r) return 0;
        assert(r->get_type() == t);
        return static_cast<_node*>(r);
    }

public:
    xsd_element* find_element_node(const string& name) const { return find_and_verify<xsd_element>(_element_map, name, xnt_element); }
    xsd_attribute* find_attribute_node(const string& name) const { return find_and_verify<xsd_attribute>(_attribute_map, name, xnt_attribute); }
    xsd_group* find_group_node(const string& name) const { return find_and_verify<xsd_group>(_group_map, name, xnt_group); }
    xsd_attribute_group* find_attribute_group_node(const string& name) const { return find_and_verify<xsd_attribute_group>(_attribute_group_map, name, xnt_attribute_group); }
    xsd_simple_type* find_simple_type_node(const string& name) const { return find_and_verify<xsd_simple_type>(_simple_type_map, name, xnt_simple_type); }
    xsd_complex_type* find_complex_type_node(const string& name) const { return find_and_verify<xsd_complex_type>(_complex_type_map, name, xnt_complex_type); }
};

struct xsd_ns_key:
    public xsd_schema
{
    xsd_ns_key(const string& str) { _namespace.assign(str); }
    xsd_ns_key(const gchar* str) { _namespace.assign(str); }
    xsd_ns_key(const gchar* str, int len) { _namespace.assign(str, len); }
};

struct xsd_loc_key:
    public xsd_schema
{
    xsd_loc_key(const string& str) { _location.assign(str); }
    xsd_loc_key(const gchar* str) { _location.assign(str); }
    xsd_loc_key(const gchar* str, int len) { _location.assign(str, len); }
};

class __gs_novtable xsd_node abstract
{
public:
    xsd_node();
    virtual ~xsd_node() {}
    virtual xsd_node_type get_type() const = 0;

public:
    void set_name(const string& str) { _name = str; }
    void set_mangling(const string& str) { _mangling = str; }
    const string& get_name() const { return _name; }
    const string& get_mangling() const { return _mangling; }
    xsd_wrapper* get_source() const { return _source; }
    void set_source(xsd_wrapper* s) { _source = s; }
    void set_translated(bool b) { _translated = b; }
    bool is_translated() const { return _translated; }

protected:
    string                  _name;
    string                  _mangling;
    xsd_wrapper*            _source;
    bool                    _translated;
};

struct xsd_node_key:
    public xsd_node
{
    xsd_node_key(const string& str) { _name.assign(str); }
    xsd_node_key(const gchar* str) { _name.assign(str); }
    xsd_node_key(const gchar* str, int len) { _name.assign(str, len); }
    xsd_node_type get_type() const override { return xnt_key; }
};

class xsd_element:
    public xsd_node
{
public:
    xsd_element();
    ~xsd_element();
    xsd_node_type get_type() const override { return xnt_element; }
    void set_ref(xsd_element* r) { _ref = r; }
    xsd_element* get_ref() const { return _ref; }
    void set_local_type(xsd_node* n) { _local_type = n; }
    xsd_node* get_local_type() const { return _local_type; }
    string& get_type_name() { return _elem_type_name; }
    string& get_entity_name() { return _elem_entity_name; }
    string& get_verificator() { return _elem_verificator; }
    string& get_assigner() { return _elem_assigner; }
    const string& const_type_name() const { return _elem_type_name; }
    const string& const_entity_name() const { return _elem_entity_name; }
    const string& const_verificator() const { return _elem_verificator; }
    const string& const_assigner() const { return _elem_assigner; }

protected:
    xsd_element*            _ref;
    xsd_node*               _local_type;
    string                  _elem_type_name;
    string                  _elem_entity_name;
    string                  _elem_verificator;
    string                  _elem_assigner;
};

class xsd_attribute:
    public xsd_node
{
public:
    xsd_attribute();
    ~xsd_attribute();
    xsd_node_type get_type() const override { return xnt_attribute; }
    void set_ref(xsd_attribute* r) { _ref = 0; }
    xsd_attribute* get_ref() const { return _ref; }
    void set_local_type(xsd_node* n) { _local_type = 0; }
    xsd_node* get_local_type() const { return _local_type; }
    string& get_type_name() { return _attr_type_name; }
    string& get_entity_name() { return _attr_entity_name; }
    string& get_verificator() { return _attr_verificator; }
    string& get_assigner() { return _attr_assigner; }
    const string& const_type_name() const { return _attr_type_name; }
    const string& const_entity_name() const { return _attr_entity_name; }
    const string& const_verificator() const { return _attr_verificator; }
    const string& const_assigner() const { return _attr_assigner; }

protected:
    xsd_attribute*          _ref;
    xsd_node*               _local_type;
    string                  _attr_type_name;
    string                  _attr_entity_name;
    string                  _attr_verificator;
    string                  _attr_assigner;
};

class xsd_group:
    public xsd_node
{
public:
    xsd_node_type get_type() const override { return xnt_group; }

protected:
};

class xsd_attribute_group:
    public xsd_node
{
public:
    xsd_node_type get_type() const override { return xnt_attribute_group; }
};

enum xsd_sts_type
{
    xst_restriction,
    xst_union,
    xst_list,
};

class __gs_novtable xsd_st_substitution abstract
{
public:
    virtual ~xsd_st_substitution() {}
    virtual xsd_sts_type get_type() const = 0;
};

class xsd_st_restriction:
    public xsd_st_substitution
{
public:
    xsd_sts_type get_type() const override { return xst_restriction; }

protected:

};

class xsd_st_union:
    public xsd_st_substitution
{
public:
    xsd_sts_type get_type() const override { return xst_union; }
};

class xsd_st_list:
    public xsd_st_substitution
{
public:
    xsd_sts_type get_type() const override { return xst_list; }
};

class xsd_simple_type:
    public xsd_node
{
public:
    xsd_simple_type();
    ~xsd_simple_type();
    xsd_node_type get_type() const override { return xnt_simple_type; }
    void set_ref(xsd_simple_type* r) { _ref = r; }
    xsd_simple_type* get_ref() const { return _ref; }
    void set_substitution(xsd_st_substitution* sub);
    bool build_assigner(string& code, const string& arg, const string& indent);

protected:
    xsd_simple_type*        _ref;
    xsd_st_substitution*    _substitution;
};

class xsd_complex_type:
    public xsd_node
{
public:
    xsd_node_type get_type() const override { return xnt_complex_type; }
    bool build_assigner(string& code, const string& arg, const string& indent) { return false; }
};

__gslib_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib