text:
/* matrix/gsl_matrix_short.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifndef __GSL_MATRIX_SHORT_H__
#define __GSL_MATRIX_SHORT_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector_short.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  short * data;
  gsl_block_short * block;
  int owner;
} gsl_matrix_short;

typedef struct
{
  gsl_matrix_short matrix;
} _gsl_matrix_short_view;

typedef _gsl_matrix_short_view gsl_matrix_short_view;

typedef struct
{
  gsl_matrix_short matrix;
} _gsl_matrix_short_const_view;

typedef const _gsl_matrix_short_const_view gsl_matrix_short_const_view;

/* Allocation */

gsl_matrix_short *
gsl_matrix_short_alloc (const size_t n1, const size_t n2);

gsl_matrix_short *
gsl_matrix_short_calloc (const size_t n1, const size_t n2);

gsl_matrix_short *
gsl_matrix_short_alloc_from_block (gsl_block_short * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

gsl_matrix_short *
gsl_matrix_short_alloc_from_matrix (gsl_matrix_short * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

gsl_vector_short *
gsl_vector_short_alloc_row_from_matrix (gsl_matrix_short * m,
                                        const size_t i);

gsl_vector_short *
gsl_vector_short_alloc_col_from_matrix (gsl_matrix_short * m,
                                        const size_t j);

void gsl_matrix_short_free (gsl_matrix_short * m);

/* Views */

_gsl_matrix_short_view
gsl_matrix_short_submatrix (gsl_matrix_short * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_short_view
gsl_matrix_short_row (gsl_matrix_short * m, const size_t i);

_gsl_vector_short_view
gsl_matrix_short_column (gsl_matrix_short * m, const size_t j);

_gsl_vector_short_view
gsl_matrix_short_diagonal (gsl_matrix_short * m);

_gsl_vector_short_view
gsl_matrix_short_subdiagonal (gsl_matrix_short * m, const size_t k);

_gsl_vector_short_view
gsl_matrix_short_superdiagonal (gsl_matrix_short * m, const size_t k);

_gsl_matrix_short_view
gsl_matrix_short_view_array (short * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_short_view
gsl_matrix_short_view_array_with_tda (short * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_short_view
gsl_matrix_short_view_vector (gsl_vector_short * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_short_view
gsl_matrix_short_view_vector_with_tda (gsl_vector_short * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_short_const_view
gsl_matrix_short_const_submatrix (const gsl_matrix_short * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_short_const_view
gsl_matrix_short_const_row (const gsl_matrix_short * m,
                            const size_t i);

_gsl_vector_short_const_view
gsl_matrix_short_const_column (const gsl_matrix_short * m,
                               const size_t j);

_gsl_vector_short_const_view
gsl_matrix_short_const_diagonal (const gsl_matrix_short * m);

_gsl_vector_short_const_view
gsl_matrix_short_const_subdiagonal (const gsl_matrix_short * m,
                                    const size_t k);

_gsl_vector_short_const_view
gsl_matrix_short_const_superdiagonal (const gsl_matrix_short * m,
                                      const size_t k);

_gsl_matrix_short_const_view
gsl_matrix_short_const_view_array (const short * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_short_const_view
gsl_matrix_short_const_view_array_with_tda (const short * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_short_const_view
gsl_matrix_short_const_view_vector (const gsl_vector_short * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_short_const_view
gsl_matrix_short_const_view_vector_with_tda (const gsl_vector_short * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

short   gsl_matrix_short_get(const gsl_matrix_short * m, const size_t i, const size_t j);
void    gsl_matrix_short_set(gsl_matrix_short * m, const size_t i, const size_t j, const short x);

short * gsl_matrix_short_ptr(gsl_matrix_short * m, const size_t i, const size_t j);
const short * gsl_matrix_short_const_ptr(const gsl_matrix_short * m, const size_t i, const size_t j);

void gsl_matrix_short_set_zero (gsl_matrix_short * m);
void gsl_matrix_short_set_identity (gsl_matrix_short * m);
void gsl_matrix_short_set_all (gsl_matrix_short * m, short x);

int gsl_matrix_short_fread (FILE * stream, gsl_matrix_short * m) ;
int gsl_matrix_short_fwrite (FILE * stream, const gsl_matrix_short * m) ;
int gsl_matrix_short_fscanf (FILE * stream, gsl_matrix_short * m);
int gsl_matrix_short_fprintf (FILE * stream, const gsl_matrix_short * m, const char * format);

int gsl_matrix_short_memcpy(gsl_matrix_short * dest, const gsl_matrix_short * src);
int gsl_matrix_short_swap(gsl_matrix_short * m1, const gsl_matrix_short * m2);

int gsl_matrix_short_swap_rows(gsl_matrix_short * m, const size_t i, const size_t j);
int gsl_matrix_short_swap_columns(gsl_matrix_short * m, const size_t i, const size_t j);
int gsl_matrix_short_swap_rowcol(gsl_matrix_short * m, const size_t i, const size_t j);
int gsl_matrix_short_transpose (gsl_matrix_short * m);
int gsl_matrix_short_transpose_memcpy (gsl_matrix_short * dest, const gsl_matrix_short * src);

short gsl_matrix_short_max (const gsl_matrix_short * m);
short gsl_matrix_short_min (const gsl_matrix_short * m);
void gsl_matrix_short_minmax (const gsl_matrix_short * m, short * min_out, short * max_out);

void gsl_matrix_short_max_index (const gsl_matrix_short * m, size_t * imax, size_t *jmax);
void gsl_matrix_short_min_index (const gsl_matrix_short * m, size_t * imin, size_t *jmin);
void gsl_matrix_short_minmax_index (const gsl_matrix_short * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

int gsl_matrix_short_isnull (const gsl_matrix_short * m);

int gsl_matrix_short_add (gsl_matrix_short * a, const gsl_matrix_short * b);
int gsl_matrix_short_sub (gsl_matrix_short * a, const gsl_matrix_short * b);
int gsl_matrix_short_mul_elements (gsl_matrix_short * a, const gsl_matrix_short * b);
int gsl_matrix_short_div_elements (gsl_matrix_short * a, const gsl_matrix_short * b);
int gsl_matrix_short_scale (gsl_matrix_short * a, const double x);
int gsl_matrix_short_add_constant (gsl_matrix_short * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_short_get_row(gsl_vector_short * v, const gsl_matrix_short * m, const size_t i);
int gsl_matrix_short_get_col(gsl_vector_short * v, const gsl_matrix_short * m, const size_t j);
int gsl_matrix_short_set_row(gsl_matrix_short * m, const size_t i, const gsl_vector_short * v);
int gsl_matrix_short_set_col(gsl_matrix_short * m, const size_t j, const gsl_vector_short * v);

extern int gsl_check_range ;

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
short
gsl_matrix_short_get(const gsl_matrix_short * m, const size_t i, const size_t j)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_short_set(gsl_matrix_short * m, const size_t i, const size_t j, const short x)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
short *
gsl_matrix_short_ptr(gsl_matrix_short * m, const size_t i, const size_t j)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (short *) (m->data + (i * m->tda + j)) ;
}

extern inline
const short *
gsl_matrix_short_const_ptr(const gsl_matrix_short * m, const size_t i, const size_t j)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const short *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_SHORT_H__ */


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_her2k (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha[2] = {0.0f, 0.1f};
   float beta = 0.1f;
   float A[] = { 0.531f, 0.721f, -0.848f, 0.826f };
   int lda = 2;
   float B[] = { -0.711f, -0.2f, -0.92f, -0.676f };
   int ldb = 2;
   float C[] = { -0.447f, 0.701f };
   int ldc = 1;
   float C_expected[] = { 0.30322f, 0.0f };
   cblas_cher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cher2k(case 1654) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cher2k(case 1654) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha[2] = {0.0f, 0.1f};
   float beta = 0.1f;
   float A[] = { 0.68f, 0.079f, 0.837f, -0.814f };
   int lda = 2;
   float B[] = { -0.986f, 0.024f, 0.584f, -0.248f };
   int ldb = 2;
   float C[] = { 0.477f, -0.551f };
   int ldc = 1;
   float C_expected[] = { 0.120103f, 0.0f };
   cblas_cher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cher2k(case 1655) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cher2k(case 1655) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha[2] = {0.0f, 0.1f};
   float beta = 0.1f;
   float A[] = { 0.354f, -0.63f, -0.85f, 0.426f };
   int lda = 1;
   float B[] = { 0.787f, -0.228f, -0.568f, 0.83f };
   int ldb = 1;
   float C[] = { 0.428f, -0.388f };
   int ldc = 1;
   float C_expected[] = { 0.0331132f, 0.0f };
   cblas_cher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cher2k(case 1656) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cher2k(case 1656) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha[2] = {0.0f, 0.1f};
   float beta = 0.1f;
   float A[] = { -0.49f, 0.224f, -0.606f, 0.46f };
   int lda = 1;
   float B[] = { -0.191f, -0.815f, 0.464f, 0.066f };
   int ldb = 1;
   float C[] = { 0.302f, 0.023f };
   int ldc = 1;
   float C_expected[] = { 0.0679396f, 0.0f };
   cblas_cher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cher2k(case 1657) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cher2k(case 1657) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 113;
   int N = 1;
   int K = 2;
   float alpha[2] = {-1.0f, 0.0f};
   float beta = 0.0f;
   float A[] = { 0.943f, 0.075f, 0.15f, -0.141f };
   int lda = 1;
   float B[] = { -0.962f, 0.422f, -0.592f, -0.789f };
   int ldb = 1;
   float C[] = { 0.728f, 0.601f };
   int ldc = 1;
   float C_expected[] = { 1.70613f, 0.0f };
   cblas_cher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cher2k(case 1658) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cher2k(case 1658) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 113;
   int N = 1;
   int K = 2;
   float alpha[2] = {-1.0f, 0.0f};
   float beta = 0.0f;
   float A[] = { -0.93f, -0.386f, 0.565f, 0.141f };
   int lda = 1;
   float B[] = { -0.801f, 0.022f, 0.558f, -0.932f };
   int ldb = 1;
   float C[] = { 0.068f, 0.501f };
   int ldc = 1;
   float C_expected[] = { -1.84059f, 0.0f };
   cblas_cher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cher2k(case 1659) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cher2k(case 1659) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 113;
   int N = 1;
   int K = 2;
   float alpha[2] = {-1.0f, 0.0f};
   float beta = 0.0f;
   float A[] = { -0.383f, 0.124f, 0.458f, -0.221f };
   int lda = 2;
   float B[] = { -0.107f, 0.199f, 0.18f, 0.122f };
   int ldb = 2;
   float C[] = { 0.896f, -0.874f };
   int ldc = 1;
   float C_expected[] = { -0.24227f, 0.0f };
   cblas_cher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cher2k(case 1660) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cher2k(case 1660) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 113;
   int N = 1;
   int K = 2;
   float alpha[2] = {-1.0f, 0.0f};
   float beta = 0.0f;
   float A[] = { 0.131f, 0.692f, 0.533f, -0.672f };
   int lda = 2;
   float B[] = { -0.435f, -0.453f, 0.195f, -0.579f };
   int ldb = 2;
   float C[] = { -0.547f, 0.736f };
   int ldc = 1;
   float C_expected[] = { -0.245124f, 0.0f };
   cblas_cher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cher2k(case 1661) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cher2k(case 1661) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta = 0.1;
   double A[] = { 0.972, -0.353, 0.712, -0.53 };
   int lda = 2;
   double B[] = { 0.787, -0.379, 0.889, 0.901 };
   int ldb = 2;
   double C[] = { 0.002, 0.266 };
   int ldc = 1;
   double C_expected[] = { -0.4278924, 0.0 };
   cblas_zher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zher2k(case 1662) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zher2k(case 1662) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta = 0.1;
   double A[] = { -0.36, 0.192, 0.539, 0.198 };
   int lda = 2;
   double B[] = { -0.673, 0.781, 0.792, 0.335 };
   int ldb = 2;
   double C[] = { 0.719, -0.339 };
   int ldc = 1;
   double C_expected[] = { -0.485009, 0.0 };
   cblas_zher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zher2k(case 1663) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zher2k(case 1663) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta = 0.1;
   double A[] = { -0.143, 0.456, 0.677, -0.522 };
   int lda = 1;
   double B[] = { 0.851, 0.196, 0.586, 0.64 };
   int ldb = 1;
   double C[] = { 0.617, 0.118 };
   int ldc = 1;
   double C_expected[] = { 0.1081226, 0.0 };
   cblas_zher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zher2k(case 1664) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zher2k(case 1664) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta = 0.1;
   double A[] = { 0.801, 0.91, 0.376, -0.006 };
   int lda = 1;
   double B[] = { -0.613, -0.758, -0.966, 0.194 };
   int ldb = 1;
   double C[] = { -0.723, -0.765 };
   int ldc = 1;
   double C_expected[] = { 0.8583678, 0.0 };
   cblas_zher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zher2k(case 1665) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zher2k(case 1665) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 113;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta = 0.1;
   double A[] = { -0.359, 0.913, 0.738, -0.227 };
   int lda = 1;
   double B[] = { 0.787, 0.745, 0.036, -0.606 };
   int ldb = 1;
   double C[] = { -0.652, -0.281 };
   int ldc = 1;
   double C_expected[] = { -0.1172608, 0.0 };
   cblas_zher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zher2k(case 1666) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zher2k(case 1666) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 113;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta = 0.1;
   double A[] = { -0.933, 0.598, 0.952, 0.25 };
   int lda = 1;
   double B[] = { -0.508, -0.461, -0.727, 0.162 };
   int ldb = 1;
   double C[] = { 0.215, 0.943 };
   int ldc = 1;
   double C_expected[] = { 0.0795166, 0.0 };
   cblas_zher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zher2k(case 1667) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zher2k(case 1667) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 113;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta = 0.1;
   double A[] = { -0.735, 0.372, -0.251, -0.168 };
   int lda = 2;
   double B[] = { 0.217, 0.863, -0.179, -0.057 };
   int ldb = 2;
   double C[] = { 0.579, -0.305 };
   int ldc = 1;
   double C_expected[] = { 0.0744312, 0.0 };
   cblas_zher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zher2k(case 1668) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zher2k(case 1668) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 113;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta = 0.1;
   double A[] = { -0.587, -0.994, -0.625, 0.681 };
   int lda = 2;
   double B[] = { -0.577, -0.014, -0.434, 0.204 };
   int ldb = 2;
   double C[] = { 0.256, 0.093 };
   int ldc = 1;
   double C_expected[] = { -0.3526202, 0.0 };
   cblas_zher2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zher2k(case 1669) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zher2k(case 1669) imag");
     };
   };
  };

}

ext: c


text:
/**
 * Licensed to the Apache Software Foundation (ASF) under one
 * or more contributor license agreements.  See the NOTICE file
 * distributed with this work for additional information
 * regarding copyright ownership.  The ASF licenses this file
 * to you under the Apache License, Version 2.0 (the
 * "License"); you may not use this file except in compliance
 * with the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#ifndef SINGA_CORE_TENSOR_TENSOR_MATH_CPP_H_
#define SINGA_CORE_TENSOR_TENSOR_MATH_CPP_H_

#include "./tensor_math.h"
//#include "./stacktrace.h"
#include <math.h>

#include <algorithm>
#include <cfloat>
#include <iostream>
#include <iterator>
#include <sstream>

#include "singa/core/common.h"
#include "singa/core/tensor.h"

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

namespace singa {

// ===================== Helper Functions =============================

// generate a traversal_info vector based on the tensor's shape for the
// traverse_next function to work
vector<int> generate_traversal_info(const Tensor &x) {
  vector<int> traversal_info = {};
  for (size_t n = 0; n < (x.shape().size() + 2); ++n) {
    traversal_info.push_back(0);
  }
  return traversal_info;
};

// generate shape multipliers
// for e.g. tensor of shape (3,3), stride (1,3) will have shape multipliers of
// (3,1)
// for e.g. tensor of shape (3,3), stride (3,1) will also have shape multipliers
// of (3,1)
// this means that the 3rd, 6th, and 9th index of the array will always be the
// starting element of their respective rows
// so we need to need use the inner stride when jumping from 1st->2nd element,
// and outer stride when jumping from 2nd->3rd
vector<int> generate_shape_multipliers(const Tensor &x) {
  Shape y_shape = x.shape();
  if (y_shape.size() == 0) {
    return {1};
  }
  vector<int> shape_multipliers = {1};
  int cumulative_product = 1;

  for (size_t n = 0; n < (y_shape.size() - 1); ++n) {
    cumulative_product = cumulative_product * y_shape[y_shape.size() - 1 - n];
    shape_multipliers.insert(shape_multipliers.begin(), cumulative_product);
  }
  return shape_multipliers;
};

// ******************************************************************************************
// CPP traversal operations (works on const declarations without modifying
// tensor variables)
// ******************************************************************************************

// this function checks whether the next index falls on a special multiplier of
// the outer shape
// so the algorithm knows when to jump over/back to a starting element of the
// outer shape
// for e.g. in [[1,4,7], [2,5,8], [3,6,9]], elements 1,2,3 are the starting
// elements of their respective rows
// this additional check only has 1 loop for 2d matrix
// but runtime performance might degrade to O(nlog(n)) for higher dimensional
// tensors
int determine_order(vector<int> &shape_multipliers, int counter) {
  for (size_t n = 0; n < (shape_multipliers.size() - 1); ++n) {
    if ((counter % shape_multipliers[n]) == 0) {
      return ((shape_multipliers.size()) - 1 - n);
    }
  }
  return 0;
};

// this function updates the base indexes with the current index after every
// single traversal step,
// can be generalized beyond 2d cases
void update_base_index(const Tensor &x, vector<int> &traversal_info) {
  for (int n = 0; n < (traversal_info[x.shape().size() + 1] + 1); ++n) {
    traversal_info[n] = traversal_info[x.shape().size()];
  }
};

// function to traverse a const strided tensor object
// it requires an additional vector, traversal_info {0,0,0,0 ...}, comprising
// (x.shape().size()+2) elements of 0
// for e.g. 2d matrix:
// index 0 and 1 store the base row and column index respectively
// index 2 stores the current index of the traversal
// index 3 stores the order of the traversal for e.g. if the order is 0,
// it means the next element can be navigated to using the innermost stride
void traverse_next(const Tensor &x, vector<int> &shape_multipliers,
                   vector<int> &traversal_info, int counter) {
  update_base_index(x, traversal_info);
  traversal_info[x.shape().size() + 1] =
      determine_order(shape_multipliers, counter);
  traversal_info[x.shape().size()] =
      traversal_info[traversal_info[x.shape().size() + 1]] +
      x.stride()[x.stride().size() - traversal_info[x.shape().size() + 1] - 1];
};

inline int next_offset(int offset, const vector<size_t> &shape,
                       const vector<int> &stride, vector<int> *index) {
  for (int k = shape.size() - 1; k >= 0; k--) {
    if (index->at(k) + 1 < int(shape.at(k))) {
      offset += stride.at(k);
      index->at(k) += 1;
      break;
    }
    index->at(k) = 0;
    offset -= stride.at(k) * (shape.at(k) - 1);
  }
  return offset;
}

template <typename DType>
void traverse_unary(const Tensor &in, Tensor *out,
                    std::function<DType(DType)> func) {
  DType *outPtr = static_cast<DType *>(out->block()->mutable_data());
  const DType *inPtr = static_cast<const DType *>(in.block()->data());
  /*
  vector<int> traversal_info = generate_traversal_info(in);
  vector<int> shape_multipliers = generate_shape_multipliers(in);

  for (size_t i = 0; i < in.Size(); i++) {
    outPtr[i] = func(inPtr[traversal_info[in.shape().size()]]);
    traverse_next(in, shape_multipliers, traversal_info, i + 1);
  }
  */
  CHECK(in.shape() == out->shape());
  if (in.stride() == out->stride()) {
    for (size_t i = 0; i < in.Size(); i++) outPtr[i] = func(inPtr[i]);
  } else {
    // LOG(INFO) << "not equal stride";
    size_t in_offset = 0, out_offset = 0;
    vector<int> in_idx(in.nDim(), 0), out_idx(out->nDim(), 0);
    for (size_t i = 0; i < Product(in.shape()); i++) {
      outPtr[out_offset] = func(inPtr[in_offset]);
      out_offset =
          next_offset(out_offset, out->shape(), out->stride(), &out_idx);
      in_offset = next_offset(in_offset, in.shape(), in.stride(), &in_idx);
    }
  }
}

template <typename DType>
void traverse_binary(const Tensor &in1, const Tensor &in2, Tensor *out,
                     std::function<DType(DType, DType)> func) {
  DType *outPtr = static_cast<DType *>(out->block()->mutable_data());
  const DType *in1Ptr = static_cast<const DType *>(in1.block()->data());
  const DType *in2Ptr = static_cast<const DType *>(in2.block()->data());
  /*
  vector<int> traversal_info_in1 = generate_traversal_info(in1);
  vector<int> traversal_info_in2 = generate_traversal_info(in2);
  vector<int> shape_multipliers_in1 = generate_shape_multipliers(in1);
  vector<int> shape_multipliers_in2 = generate_shape_multipliers(in2);

  for (size_t i = 0; i < in1.Size(); i++) {
    outPtr[i] = func(in1Ptr[traversal_info_in1[in1.shape().size()]],
                     in2Ptr[traversal_info_in2[in2.shape().size()]]);
    traverse_next(in1, shape_multipliers_in1, traversal_info_in1, i + 1);
    traverse_next(in2, shape_multipliers_in2, traversal_info_in2, i + 1);
  }
  */
  auto prod = Product(in1.shape());
  CHECK(in1.shape() == out->shape());
  CHECK(in2.shape() == out->shape());
  if ((in1.stride() == out->stride()) && (in2.stride() == in1.stride())) {
    for (size_t i = 0; i < prod; i++) outPtr[i] = func(in1Ptr[i], in2Ptr[i]);
  } else {
    /*
    LOG(INFO) << "not equal stride";
    std::ostringstream s1, s2, s3, s4, s5, s6;
    std::copy(in1.stride().begin(), in1.stride().end(),
    std::ostream_iterator<int>(s1, ", "));
    std::copy(in2.stride().begin(), in2.stride().end(),
    std::ostream_iterator<int>(s2, ", "));
    std::copy(out->stride().begin(), out->stride().end(),
    std::ostream_iterator<int>(s3, ", "));

    std::copy(in1.shape().begin(), in1.shape().end(),
    std::ostream_iterator<int>(s4, ", "));
    std::copy(in2.shape().begin(), in2.shape().end(),
    std::ostream_iterator<int>(s5, ", "));
    std::copy(out->shape().begin(), out->shape().end(),
    std::ostream_iterator<int>(s6, ", "));

    LOG(INFO) << s1.str() << ": " << s4.str();
    LOG(INFO) << s2.str() << ": " << s5.str();
    LOG(INFO) << s3.str() << ": " << s6.str();
    LOG(INFO) << Backtrace();
    */

    size_t in1_offset = 0, in2_offset = 0, out_offset = 0;
    vector<int> in1_idx(in1.nDim(), 0), in2_idx(in2.nDim(), 0),
        out_idx(out->nDim(), 0);
    for (size_t i = 0; i < prod; i++) {
      outPtr[out_offset] = func(in1Ptr[in1_offset], in2Ptr[in2_offset]);
      out_offset =
          next_offset(out_offset, out->shape(), out->stride(), &out_idx);
      in1_offset = next_offset(in1_offset, in1.shape(), in1.stride(), &in1_idx);
      in2_offset = next_offset(in2_offset, in2.shape(), in2.stride(), &in2_idx);
      // LOG(INFO) <<  in1_offset << ", " << in2_offset << ", " << out_offset;
    }
  }
}

// ******************************************************************************************
// traversal operations end
// ******************************************************************************************

// ===================== CUDA Functions =============================

template <>
void Abs<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  traverse_unary<float>(in, out, [](float x) { return fabs(x); });
}

template <>
void Erf<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  traverse_unary<float>(in, out, [](float x) { return erff(x); });
}

template <>
void CastCopy<float, half_float::half, lang::Cpp>(const Tensor *src,
                                                  Tensor *dst, Context *ctx) {
  half_float::half *dst_array =
      static_cast<half_float::half *>(dst->block()->mutable_data());
  const float *src_array = static_cast<const float *>(src->block()->data());
  for (int i = 0; i < dst->Size(); ++i)
    dst_array[i] = static_cast<half_float::half>(src_array[i]);
}

template <>
void CastCopy<half_float::half, float, lang::Cpp>(const Tensor *src,
                                                  Tensor *dst, Context *ctx) {
  float *dst_array = static_cast<float *>(dst->block()->mutable_data());
  const half_float::half *src_array =
      static_cast<const half_float::half *>(src->block()->data());
  for (int i = 0; i < dst->Size(); ++i)
    dst_array[i] = static_cast<float>(src_array[i]);
}

template <>
void CastCopy<float, int, lang::Cpp>(const Tensor *src, Tensor *dst,
                                     Context *ctx) {
  int *dst_array = static_cast<int *>(dst->block()->mutable_data());
  const float *src_array = static_cast<const float *>(src->block()->data());
  for (int i = 0; i < dst->Size(); ++i) dst_array[i] = (int)src_array[i];
}

template <>
void CastCopy<int, float, lang::Cpp>(const Tensor *src, Tensor *dst,
                                     Context *ctx) {
  float *dst_array = static_cast<float *>(dst->block()->mutable_data());
  const int *src_array = static_cast<const int *>(src->block()->data());
  for (int i = 0; i < dst->Size(); ++i) dst_array[i] = (float)src_array[i];
}

template <>
void Ceil<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  traverse_unary<float>(in, out, [](float x) { return std::ceil(x); });
}

template <>
void Floor<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  traverse_unary<float>(in, out, [](float x) { return std::floor(x); });
}

template <>
void Round<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  traverse_unary<float>(in, out, [](float x) { return std::round(x); });
}

template <>
void RoundE<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  traverse_unary<float>(in, out, [](float x) {
    float doub = x * 2;
    if (ceilf(doub) == doub) {
      return std::round(x / 2) * 2;
    } else {
      return std::round(x);
    }
  });
}

#ifdef USE_DNNL
template <>
void SoftMax<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto md = dnnl::memory::desc({static_cast<long long>(in.shape()[0]),
                                static_cast<long long>(in.shape()[1])},
                               dnnl::memory::data_type::f32,
                               dnnl::memory::format_tag::ab);
  auto in_mem = dnnl::memory(md, ctx->dnnl_engine, in.block()->mutable_data());
  auto out_mem =
      dnnl::memory(md, ctx->dnnl_engine, out->block()->mutable_data());

  auto softmax_desc =
      dnnl::softmax_forward::desc(dnnl::prop_kind::forward_scoring, md, 1);
  auto softmax_prim_desc =
      dnnl::softmax_forward::primitive_desc(softmax_desc, ctx->dnnl_engine);
  auto softmax = dnnl::softmax_forward(softmax_prim_desc);
  softmax.execute(ctx->dnnl_stream,
                  {{DNNL_ARG_SRC, in_mem}, {DNNL_ARG_DST, out_mem}});
  ctx->dnnl_stream.wait();
}

template <>
void SoftMaxBackward<float, lang::Cpp>(const Tensor &in, Tensor *out,
                                       const Tensor &fdout, Context *ctx) {
  auto md = dnnl::memory::desc({static_cast<long long>(in.shape()[0]),
                                static_cast<long long>(in.shape()[1])},
                               dnnl::memory::data_type::f32,
                               dnnl::memory::format_tag::ab);
  auto in_mem = dnnl::memory(md, ctx->dnnl_engine, in.block()->mutable_data());
  auto fdout_mem =
      dnnl::memory(md, ctx->dnnl_engine, fdout.block()->mutable_data());
  auto out_mem =
      dnnl::memory(md, ctx->dnnl_engine, out->block()->mutable_data());

  auto softmax_desc =
      dnnl::softmax_forward::desc(dnnl::prop_kind::forward_scoring, md, 1);
  auto softmax_prim_desc =
      dnnl::softmax_forward::primitive_desc(softmax_desc, ctx->dnnl_engine);

  auto softmaxbwd_desc = dnnl::softmax_backward::desc(md, md, 1);
  auto softmaxbwd_prim_desc = dnnl::softmax_backward::primitive_desc(
      softmaxbwd_desc, ctx->dnnl_engine, softmax_prim_desc);
  auto softmaxbwd = dnnl::softmax_backward(softmaxbwd_prim_desc);
  softmaxbwd.execute(ctx->dnnl_stream, {{DNNL_ARG_DIFF_SRC, out_mem},
                                        {DNNL_ARG_DIFF_DST, in_mem},
                                        {DNNL_ARG_DST, fdout_mem}});
  ctx->dnnl_stream.wait();
}
#else
// native Softmax without DNNL
template <>
void SoftMax<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  CHECK_LE(in.nDim(), 2u)
      << "Axis is required for SoftMax on multi dimemsional tensor";
  out->CopyData(in);
  size_t nrow = 1, ncol = in.Size(), size = ncol;
  if (in.nDim() == 2u) {
    nrow = in.shape(0);
    ncol = size / nrow;
    out->Reshape(Shape{nrow, ncol});
  }
  Tensor tmp = RowMax(*out);
  SubColumn(tmp, out);
  Exp(*out, out);

  SumColumns(*out, &tmp);
  DivColumn(tmp, out);
  out->Reshape(in.shape());
}
#endif  // USE_DNNL

template <>
void Add<float, lang::Cpp>(const Tensor &in, const float x, Tensor *out,
                           Context *ctx) {
  auto add_lambda = [&x](float a) { return (a + x); };
  traverse_unary<float>(in, out, add_lambda);
}

template <>
void Add<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                           Context *ctx) {
  // CHECK_EQ(ctx->stream, nullptr);
  auto add_lambda_binary = [](float a, float b) { return (a + b); };
  traverse_binary<float>(in1, in2, out, add_lambda_binary);
}

template <>
void Clamp<float, lang::Cpp>(const float low, const float high,
                             const Tensor &in, Tensor *out, Context *ctx) {
  auto clamp_lambda = [&low, &high](float a) {
    if (a < low) {
      return low;
    } else if (a > high) {
      return high;
    } else {
      return a;
    }
  };
  traverse_unary<float>(in, out, clamp_lambda);
}

template <>
void Div<float, lang::Cpp>(const float x, const Tensor &in, Tensor *out,
                           Context *ctx) {
  auto const_div = [&x](float a) {
    CHECK_NE(a, 0.f);
    return x / a;
  };
  traverse_unary<float>(in, out, const_div);
}

template <>
void Div<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                           Context *ctx) {
  auto binary_div = [](float a, float b) {
    CHECK_NE(b, 0.f);
    return a / b;
  };
  traverse_binary<float>(in1, in2, out, binary_div);
}

template <>
void EltwiseMult<float, lang::Cpp>(const Tensor &in, const float x, Tensor *out,
                                   Context *ctx) {
  auto eltwisemult_lambda = [&x](float a) { return (a * x); };
  traverse_unary<float>(in, out, eltwisemult_lambda);
}

template <>
void EltwiseMult<float, lang::Cpp>(const Tensor &in1, const Tensor &in2,
                                   Tensor *out, Context *ctx) {
  auto eltwisemult_lambda_binary = [](float a, float b) { return (a * b); };
  traverse_binary<float>(in1, in2, out, eltwisemult_lambda_binary);
}

template <>
void ReLUBackward<float, lang::Cpp>(const Tensor &in1, const Tensor &in2,
                                    Tensor *out, Context *ctx) {
  auto relubackward_lambda = [](float a, float b) { return (b > 0) ? a : 0.f; };
  traverse_binary<float>(in1, in2, out, relubackward_lambda);
}

template <>
void Exp<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  traverse_unary<float>(in, out, [](float x) { return exp(x); });
}

template <>
void GE<float, lang::Cpp>(const Tensor &in, const float x, Tensor *out,
                          Context *ctx) {
  auto ge_lambda = [&x](float a) { return (a >= x) ? 1.f : 0.f; };
  traverse_unary<float>(in, out, ge_lambda);
}

template <>
void GE<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                          Context *ctx) {
  auto ge_lambda_binary = [](float a, float b) { return (a >= b) ? 1.f : 0.f; };
  traverse_binary<float>(in1, in2, out, ge_lambda_binary);
}

template <>
void GE<int, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                        Context *ctx) {
  auto ge_lambda_binary = [](int a, int b) { return (a >= b) ? 1.f : 0.f; };
  traverse_binary<int>(in1, in2, out, ge_lambda_binary);
}

template <>
void GT<float, lang::Cpp>(const Tensor &in, const float x, Tensor *out,
                          Context *ctx) {
  auto gt_lambda = [&x](float a) { return (a > x) ? 1.f : 0.f; };
  traverse_unary<float>(in, out, gt_lambda);
}

template <>
void GT<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                          Context *ctx) {
  auto gt_lambda_binary = [](float a, float b) { return (a > b) ? 1.f : 0.f; };
  traverse_binary<float>(in1, in2, out, gt_lambda_binary);
}

template <>
void GT<int, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                        Context *ctx) {
  auto gt_lambda_binary = [](int a, int b) { return (a > b) ? 1.f : 0.f; };
  traverse_binary<int>(in1, in2, out, gt_lambda_binary);
}

template <>
void LE<float, lang::Cpp>(const Tensor &in, const float x, Tensor *out,
                          Context *ctx) {
  auto le_lambda = [&x](float a) { return (a <= x) ? 1.f : 0.f; };
  traverse_unary<float>(in, out, le_lambda);
}

template <>
void LE<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                          Context *ctx) {
  auto le_lambda_binary = [](float a, float b) { return (a <= b) ? 1.f : 0.f; };
  traverse_binary<float>(in1, in2, out, le_lambda_binary);
}

template <>
void LE<int, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                        Context *ctx) {
  auto le_lambda_binary = [](int a, int b) { return (a <= b) ? 1.f : 0.f; };
  traverse_binary<int>(in1, in2, out, le_lambda_binary);
}

template <>
void Log<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto ulog = [](float a) {
    CHECK_GT(a, 0.f);
    return log(a);
  };
  traverse_unary<float>(in, out, ulog);
}

template <>
void LT<float, lang::Cpp>(const Tensor &in, const float x, Tensor *out,
                          Context *ctx) {
  auto lt_lambda = [&x](float a) { return (a < x) ? 1.f : 0.f; };
  traverse_unary<float>(in, out, lt_lambda);
}

template <>
void LT<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                          Context *ctx) {
  auto lt_lambda_binary = [](float a, float b) { return (a < b) ? 1.f : 0.f; };
  traverse_binary<float>(in1, in2, out, lt_lambda_binary);
}

template <>
void LT<int, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                        Context *ctx) {
  auto lt_lambda_binary = [](int a, int b) { return (a < b) ? 1.f : 0.f; };
  traverse_binary<int>(in1, in2, out, lt_lambda_binary);
}

template <>
void EQ<float, lang::Cpp>(const Tensor &in, const float x, Tensor *out,
                          Context *ctx) {
  auto eq_lambda = [&x](float a) { return (a == x) ? 1.f : 0.f; };
  traverse_unary<float>(in, out, eq_lambda);
}

template <>
void EQ<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                          Context *ctx) {
  auto eq_lambda_binary = [](float a, float b) { return (a == b) ? 1.f : 0.f; };
  traverse_binary<float>(in1, in2, out, eq_lambda_binary);
}

template <>
void EQ<int, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                        Context *ctx) {
  auto eq_lambda_binary = [](int a, int b) { return (a == b) ? 1.f : 0.f; };
  traverse_binary<int>(in1, in2, out, eq_lambda_binary);
}

template <>
void Pow<float, lang::Cpp>(const Tensor &in, const float x, Tensor *out,
                           Context *ctx) {
  traverse_unary<float>(in, out, [x](float y) { return pow(y, x); });
}

template <>
void Pow<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                           Context *ctx) {
  auto pow_lambda_binary = [](float a, float b) { return pow(a, b); };
  traverse_binary<float>(in1, in2, out, pow_lambda_binary);
}

template <>
void ReLU<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto relu_lambda = [](float a) { return (a >= 0.f) ? a : 0.f; };
  traverse_unary<float>(in, out, relu_lambda);
}

template <>
void Set<float, lang::Cpp>(const float x, Tensor *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  for (size_t i = 0; i < out->Size(); i++) outPtr[i] = x;
}

template <>
void Set<int, lang::Cpp>(const int x, Tensor *out, Context *ctx) {
  int *outPtr = static_cast<int *>(out->block()->mutable_data());
  for (size_t i = 0; i < out->Size(); i++) outPtr[i] = x;
}

template <>
void Set<half_float::half, lang::Cpp>(const half_float::half x, Tensor *out,
                                      Context *ctx) {
  half_float::half *outPtr =
      static_cast<half_float::half *>(out->block()->mutable_data());
  for (size_t i = 0; i < out->Size(); i++) outPtr[i] = x;
}

template <>
void Sigmoid<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto sigmoid_lambda = [](float a) { return 1.f / (1.f + exp(-a)); };
  traverse_unary<float>(in, out, sigmoid_lambda);
}

template <>
void Sign<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto sign_lambda = [](float a) { return (a > 0) - (a < 0); };
  traverse_unary<float>(in, out, sign_lambda);
}

template <>
void SoftPlus<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto softplus_lambda = [](float a) { return log(1.f + exp(a)); };
  traverse_unary<float>(in, out, softplus_lambda);
}

template <>
void SoftSign<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto softsign_lambda = [](float a) { return a / (1.f + fabs(a)); };
  traverse_unary<float>(in, out, softsign_lambda);
}

template <>
void Sqrt<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto usqrt = [](float a) {
    CHECK_GE(a, 0.f);
    return sqrt(a);
  };
  traverse_unary<float>(in, out, usqrt);
}

template <>
void Sub<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                           Context *ctx) {
  // CHECK_EQ(ctx->stream, nullptr);
  auto sub_lambda_binary = [](float a, float b) { return (a - b); };
  traverse_binary<float>(in1, in2, out, sub_lambda_binary);
}

// sum all elements of input into out
// TODO(wangwei) optimize using omp
template <>
void Sum<float, lang::Cpp>(const Tensor &in, float *out, Context *ctx) {
  float s = 0.f;
  const float *inPtr = static_cast<const float *>(in.block()->data());
  for (size_t i = 0; i < in.Size(); i++) {
    s += inPtr[i];
  }
  *out = s;
}

#define GenUnaryTensorCppFn(fn, cppfn)                                     \
  template <>                                                              \
  void fn<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) { \
    auto fn_lambda = [](float a) { return cppfn(a); };                     \
    traverse_unary<float>(in, out, fn_lambda);                             \
  }

GenUnaryTensorCppFn(Cos, cos);
GenUnaryTensorCppFn(Cosh, cosh);
GenUnaryTensorCppFn(Acos, acos);
GenUnaryTensorCppFn(Acosh, acosh);
GenUnaryTensorCppFn(Sin, sin);
GenUnaryTensorCppFn(Sinh, sinh);
GenUnaryTensorCppFn(Asin, asin);
GenUnaryTensorCppFn(Asinh, asinh);
GenUnaryTensorCppFn(Tan, tan);
GenUnaryTensorCppFn(Tanh, tanh);
GenUnaryTensorCppFn(Atan, atan);
GenUnaryTensorCppFn(Atanh, atanh);

template <>
void Transform<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto identity = [](float a) { return a; };
  traverse_unary<float>(in, out, identity);
}

template <>
void Transform<int, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  auto identity = [](int a) { return a; };
  traverse_unary<int>(in, out, identity);
}

template <>
void Transform<half_float::half, lang::Cpp>(const Tensor &in, Tensor *out,
                                            Context *ctx) {
  auto identity = [](half_float::half a) { return a; };
  traverse_unary<half_float::half>(in, out, identity);
}

template <>
void Bernoulli<float, lang::Cpp>(const float p, Tensor *out, Context *ctx) {
  std::bernoulli_distribution distribution(p);
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  for (size_t i = 0; i < out->Size(); i++) {
    outPtr[i] = distribution(ctx->random_generator) ? 1.0f : 0.0f;
  }
}

template <>
void Gaussian<float, lang::Cpp>(const float mean, const float std, Tensor *out,
                                Context *ctx) {
  std::normal_distribution<float> distribution(mean, std);
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  for (size_t i = 0; i < out->Size(); i++) {
    outPtr[i] = static_cast<float>(distribution(ctx->random_generator));
  }
}

template <>
void Gaussian<half_float::half, lang::Cpp>(const half_float::half mean,
                                           const half_float::half std,
                                           Tensor *out, Context *ctx) {
  Tensor tmp(out->shape(), out->device(), kFloat32);
  Gaussian<float, lang::Cpp>(static_cast<float>(mean), static_cast<float>(std),
                             &tmp, ctx);
  CastCopy<float, half_float::half, lang::Cpp>(&tmp, out, ctx);
}

template <>
void Uniform<float, lang::Cpp>(const float low, const float high, Tensor *out,
                               Context *ctx) {
  std::uniform_real_distribution<float> distribution(low, high);
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  for (size_t i = 0; i < out->Size(); i++) {
    outPtr[i] = static_cast<float>(distribution(ctx->random_generator));
  }
}

// ====================Blas operations======================================

// warning, this function has block M overwritting to block M itself
template <>
void DGMM<float, lang::Cpp>(const bool side_right, const Tensor &M,
                            const Tensor &v, Tensor *out, Context *ctx) {
  const float *MPtr = static_cast<const float *>(M.block()->data());
  const float *vPtr = static_cast<const float *>(v.block()->data());
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  const size_t nrow = M.shape(0);
  const size_t ncol = M.shape(1);

  if (side_right) {
    for (size_t r = 0; r < nrow; r++) {
      size_t in_offset = M.stride()[0] * r, out_offset = out->stride()[0] * r;
      for (size_t c = 0; c < ncol; c++) {
        outPtr[out_offset] = MPtr[in_offset] * vPtr[c];
        in_offset += M.stride()[1];
        out_offset += out->stride()[1];
      }
    }
  } else {
    for (size_t r = 0; r < nrow; r++) {
      size_t in_offset = M.stride()[0] * r, out_offset = out->stride()[0] * r;
      for (size_t c = 0; c < ncol; c++) {
        outPtr[out_offset] = MPtr[in_offset] * vPtr[r];
        in_offset += M.stride()[1];
        out_offset += out->stride()[1];
      }
    }
  }
}

#ifdef USE_CBLAS
template <>
void Amax<float, lang::Cpp>(const Tensor &in, size_t *out, Context *ctx) {
  const float *inPtr = static_cast<const float *>(in.block()->data());
  *out = cblas_isamax(in.Size(), inPtr, 1);  // not using strided traversal
}

template <>
void Asum<float, lang::Cpp>(const Tensor &in, float *out, Context *ctx) {
  const float *inPtr = static_cast<const float *>(in.block()->data());
  *out = cblas_sasum(in.Size(), inPtr, 1);  // not using strided traversal
}

// template <>
// void Axpy<float, lang::Cpp>(const float alpha,
//                             const Tensor& in, Tensor *out, Context *ctx) {
//   //check input tensor for strides first
//   if (in.stride() == out->stride()) {
//     const float *inPtr = static_cast<const float *>(in.block()->data());
//     float *outPtr = static_cast<float *>(out->block()->mutable_data());
//     cblas_saxpy(in.Size(), alpha, inPtr, 1, outPtr, 1);
//   } else {
//     //LOG(FATAL) << "Axpy, input and output strides do not match." ;
//     EltwiseMult<float, lang::Cpp>(in, alpha, out, ctx);
//   }
// }

template <>
void Axpy<float, lang::Cpp>(const float alpha, const Tensor &in, Tensor *out,
                            Context *ctx) {
  // check input tensor for strides first
  const float *inPtr = static_cast<const float *>(in.block()->data());
  float *outPtr = static_cast<float *>(out->block()->mutable_data());

  if (in.stride() == out->stride()) {
    cblas_saxpy(in.Size(), alpha, inPtr, 1, outPtr, 1);
  } else {
    // LOG(FATAL) << "Axpy, input and output strides do not match." ;
    Tensor t(in.shape(), in.device(), in.data_type());
    EltwiseMult<float, lang::Cpp>(in, alpha, &t, ctx);
    float *tPtr = static_cast<float *>(t.block()->mutable_data());
    cblas_saxpy(in.Size(), 1, tPtr, 1, outPtr, 1);
  }
}

// template <>
// void Axpy<float, lang::Cpp>(const float alpha,
//                            const Tensor& in, Tensor *out, Context *ctx) {
//  //check input tensor for strides first
//  if (in.stride() == out->stride()) {
//    const float *inPtr = static_cast<const float *>(in.block()->data());
//    float *outPtr = static_cast<float *>(out->block()->mutable_data());
//    cblas_saxpy(in.Size(), alpha, inPtr, 1, outPtr, 1);
//  } else if(out->transpose()) {
//    LOG(FATAL) << "output is already transposed." ;
//  } else {
//    LOG(FATAL) << "Axpy, input and output strides do not match." ;
//  }
// }

template <>
void Dot<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, float *out,
                           Context *ctx) {
  // check input tensor for strides first
  if (!(in1.transpose()) && !(in2.transpose())) {
    const float *in1Ptr = static_cast<const float *>(in1.block()->data());
    const float *in2Ptr = static_cast<const float *>(in2.block()->data());
    *out = cblas_sdot(in1.Size(), in1Ptr, 1, in2Ptr, 1);
  } else {
    LOG(FATAL) << "Dot, one of the input is tranposed. Not implemented yet.";
  }
}
template <>
void Dot<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, Tensor *out,
                           Context *ctx) {
  // check input tensor for strides first
  if (!(in1.transpose()) && !(in2.transpose())) {
    const float *in1Ptr = static_cast<const float *>(in1.block()->data());
    const float *in2Ptr = static_cast<const float *>(in2.block()->data());
    float *outPtr = static_cast<float *>(out->block()->mutable_data());
    *outPtr = cblas_sdot(in1.Size(), in1Ptr, 1, in2Ptr, 1);
  } else {
    LOG(FATAL) << "Dot, one of the input is tranposed. Not implemented yet.";
  }
}

template <>
void Scale<float, lang::Cpp>(const float x, Tensor *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  cblas_sscal(out->Size(), x, outPtr, 1);  // not using strided traversal
}

template <>
void Nrm2<float, lang::Cpp>(const Tensor &in, float *out, Context *ctx) {
  const float *inPtr = static_cast<const float *>(in.block()->data());
  *out = cblas_snrm2(in.Size(), inPtr, 1);  // not using strided traversal
}

template <>
void GEMV<float, lang::Cpp>(const float alpha, const Tensor &A, const Tensor &v,
                            const float beta, Tensor *out, Context *ctx) {
  const float *APtr = static_cast<const float *>(A.block()->data());
  const float *vPtr = static_cast<const float *>(v.block()->data());
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  const size_t m = A.shape()[0];
  const size_t n = A.shape()[1];
  if (A.transpose()) {
    cblas_sgemv(CblasRowMajor, CblasTrans, n, m, alpha, APtr, m, vPtr, 1, beta,
                outPtr, 1);
  } else {
    cblas_sgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, APtr, n, vPtr, 1,
                beta, outPtr, 1);
  }
}

template <>
void GEMM<float, lang::Cpp>(const float alpha, const Tensor &A, const Tensor &B,
                            const float beta, Tensor *C, Context *ctx) {
  auto transA = A.transpose();
  auto transa = transA ? CblasTrans : CblasNoTrans;
  auto transB = B.transpose();
  auto transb = transB ? CblasTrans : CblasNoTrans;
  const size_t nrowA = A.shape()[0];
  const size_t ncolA = A.shape()[1];
  const size_t ncolB = B.shape()[1];
  auto lda = transA ? nrowA : ncolA;
  auto ldb = transB ? ncolA : ncolB;
  auto ldc = ncolB;
  const float *APtr = static_cast<const float *>(A.block()->data());
  const float *BPtr = static_cast<const float *>(B.block()->data());
  float *CPtr = static_cast<float *>(C->block()->mutable_data());
  cblas_sgemm(CblasRowMajor, transa, transb, nrowA, ncolB, ncolA, alpha, APtr,
              lda, BPtr, ldb, beta, CPtr, ldc);
}

/*
 * implement matmul for 3d 4d tensor
 *   simulate cblas_sgemm_batch();
 *   which is only available in intel cblas
 */
template <>
void GEMMBatched<float, lang::Cpp>(const float alpha, const Tensor &A,
                                   const Tensor &B, const float beta, Tensor *C,
                                   Context *ctx) {
  const float *APtr = static_cast<const float *>(A.block()->data());
  const float *BPtr = static_cast<const float *>(B.block()->data());
  float *CPtr = static_cast<float *>(C->block()->mutable_data());

  auto transA = A.transpose();
  auto transa = transA ? CblasTrans : CblasNoTrans;
  auto transB = B.transpose();
  auto transb = transB ? CblasTrans : CblasNoTrans;

  const size_t ncolB = B.shape().end()[-1];
  const size_t nrowA = A.shape().end()[-2];
  const size_t ncolA = A.shape().end()[-1];

  auto lda = transA ? nrowA : ncolA;
  auto ldb = transB ? ncolA : ncolB;
  auto ldc = ncolB;
  const int group_count = 1;

  size_t group_size = A.shape()[0];                // 3d
  if (A.nDim() == 4u) group_size *= A.shape()[1];  // 4d

  auto matrix_stride_A = A.shape().end()[-1] * A.shape().end()[-2];
  auto matrix_stride_B = B.shape().end()[-1] * B.shape().end()[-2];
  auto matrix_stride_C = C->shape().end()[-1] * C->shape().end()[-2];
  auto offset_A = 0;
  auto offset_B = 0;
  auto offset_C = 0;

  for (int i = 0; i < group_size; i++) {
    cblas_sgemm(CblasRowMajor, transa, transb, nrowA, ncolB, ncolA, alpha,
                APtr + offset_A, lda, BPtr + offset_B, ldb, beta,
                CPtr + offset_C, ldc);
    offset_A += matrix_stride_A;
    offset_B += matrix_stride_B;
    offset_C += matrix_stride_C;
  }
}

#else

template <>
void Amax<float, lang::Cpp>(const Tensor &in, size_t *out, Context *ctx) {
  size_t maxPos = 0;
  float maxVal = 0;
  const float *inPtr = static_cast<const float *>(in.block()->data());
  for (size_t i = 0; i < in.Size(); i++) {  // not using strided traversal
    if (i == 0) {
      maxVal = inPtr[i];
    } else if (inPtr[i] > maxVal) {
      maxVal = inPtr[i];
      maxPos = i;
    }
  }
  *out = maxPos;
}

template <>
void Amin<float, lang::Cpp>(const Tensor &in, size_t *out, Context *ctx) {
  size_t minPos = 0;
  float minVal = 0;
  const float *inPtr = static_cast<const float *>(in.block()->data());
  for (size_t i = 0; i < in.Size(); i++) {  // not using strided traversal
    if (i == 0) {
      minVal = inPtr[i];
    } else if (inPtr[i] > minVal) {
      minVal = inPtr[i];
      minPos = i;
    }
  }
  *out = minPos;
}

template <>
void Asum<float, lang::Cpp>(const Tensor &in, float *out, Context *ctx) {
  float sum = 0;
  const float *inPtr = static_cast<const float *>(in.block()->data());
  for (size_t i = 0; i < in.Size(); i++) {
    sum += fabs(inPtr[i]);  // not using strided traversal
  }
}

template <>
void Axpy<float, lang::Cpp>(const float alpha, const Tensor &in, Tensor *out,
                            Context *ctx) {
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  const float *inPtr = static_cast<const float *>(in.block()->data());
  vector<int> traversal_info = generate_traversal_info(in);
  vector<int> shape_multipliers = generate_shape_multipliers(in);

  for (size_t i = 0; i < in.Size(); i++) {
    outPtr[i] += alpha * inPtr[traversal_info[in.shape().size()]];
    traverse_next(in, shape_multipliers, traversal_info, i + 1);
  }
}

template <>
void Scale<float, lang::Cpp>(const float x, Tensor *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  for (size_t i = 0; i < out->Size(); i++) {
    outPtr[i] *= x;  // not using strided traversal
  }
}

template <>
void Dot<float, lang::Cpp>(const Tensor &in1, const Tensor &in2, float *out,
                           Context *ctx) {
  float sum = 0;
  // const float *in1Ptr = static_cast<const float *>(in1.data());
  // const float *in2Ptr = static_cast<const float *>(in2.data());
  // for (size_t i = 0; i < in.Size(); i++) {
  //   sum += in1Ptr[i] * in2Ptr[i];
  // }
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  const float *in1Ptr = static_cast<const float *>(in1.block()->data());
  const float *in2Ptr = static_cast<const float *>(in2.block()->data());
  vector<int> traversal_info_in1 = generate_traversal_info(in1);
  vector<int> traversal_info_in2 = generate_traversal_info(in2);
  vector<int> shape_multipliers_in1 = generate_shape_multipliers(in1);
  vector<int> shape_multipliers_in2 = generate_shape_multipliers(in2);

  for (size_t i = 0; i < in1.Size(); i++) {
    sum += in1Ptr[traversal_info_in1[in1.shape().size()]] *
           in2Ptr[traversal_info_in2[in2.shape().size()]];
    traverse_next(in1, shape_multipliers_in1, traversal_info_in1, i + 1);
    traverse_next(in2, shape_multipliers_in2, traversal_info_in2, i + 1);
  }
}

template <>
void GEMV<float, lang::Cpp>(const float alpha, const Tensor &A, const Tensor &v,
                            const float beta, Tensor *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  const float *APtr = static_cast<const float *>(A.block()->data());
  const float *vPtr = static_cast<const float *>(v.block()->data());
  bool trans = A.transpose();
  const size_t m = A.shape(0);
  const size_t n = A.shape(1);
  for (size_t r = 0; r < m; r++) {
    float sum = 0;
    for (size_t c = 0; c < n; c++) {
      size_t idx = trans ? c * m + r : r * n + c;
      sum += APtr[idx] * vPtr[c];
    }
    outPtr[r] = alpha * sum + beta * outPtr[r];
  }
}

#endif  // USE_CBLAS
template <>
void ComputeCrossEntropy<float, lang::Cpp>(bool int_target,
                                           const size_t batchsize,
                                           const size_t dim, const Tensor &p,
                                           const Tensor &t, Tensor *loss,
                                           Context *ctx) {
  const float *pPtr = static_cast<const float *>(p.block()->data());
  const int *tPtr = static_cast<const int *>(t.block()->data());
  float *lossPtr = static_cast<float *>(loss->block()->mutable_data());
  if (int_target) {
    for (size_t i = 0; i < batchsize; i++) {
      int truth_idx = tPtr[i];
      CHECK_GE(truth_idx, 0);
      float prob_of_truth = pPtr[i * dim + truth_idx];
      lossPtr[i] = -std::log((std::max)(prob_of_truth, FLT_MIN));
    }
  } else {
    for (size_t i = 0; i < batchsize; i++) {
      float sum = 0.f;
      for (size_t j = 0; j < dim; j++) {
        sum += tPtr[i * dim + j];
      }
      float loss_value = 0.f;
      for (size_t j = 0, offset = i * dim; j < dim; j++, offset++) {
        loss_value -=
            tPtr[offset] / sum * std::log((std::max)(pPtr[offset], FLT_MIN));
      }
      lossPtr[i] = loss_value;
    }
  }
}

template <>
void SoftmaxCrossEntropyBwd<float, lang::Cpp>(bool int_target,
                                              const size_t batchsize,
                                              const size_t dim, const Tensor &p,
                                              const Tensor &t, Tensor *grad,
                                              Context *ctx) {
  CHECK_EQ(p.block(), grad->block())
      << "Use the same pointer to optimize performance";
  // const float* pPtr = static_cast<const float*>(p->data());
  const int *tPtr = static_cast<const int *>(t.block()->data());
  float *gradPtr = static_cast<float *>(grad->block()->mutable_data());

  if (int_target) {
    for (size_t i = 0; i < batchsize; i++) {
      int truth_idx = static_cast<int>(tPtr[i]);
      CHECK_GE(truth_idx, 0);
      gradPtr[i * dim + truth_idx] -= 1.0;
    }
  } else {
    for (size_t i = 0; i < batchsize; i++) {
      float sum = 0.f;
      for (size_t j = 0; j < dim; j++) {
        sum += tPtr[i * dim + j];
      }
      for (size_t j = 0, offset = i * dim; j < dim; j++, offset++) {
        gradPtr[offset] -= tPtr[offset] / sum;
      }
    }
  }
}

template <>
void RowMax<float, lang::Cpp>(const Tensor &in, Tensor *out, Context *ctx) {
  const float *inPtr = static_cast<const float *>(in.block()->data());
  float *outPtr = static_cast<float *>(out->block()->mutable_data());
  const size_t nrow = in.shape()[0];
  const size_t ncol = in.shape()[1];
  vector<int> traversal_info = generate_traversal_info(in);
  vector<int> shape_multipliers = generate_shape_multipliers(in);

  for (size_t r = 0; r < nrow; r++) {
    int counter_offset = (r * ncol);
    float maxval = 0;
    for (size_t c = 0; c < ncol; c++) {
      maxval = (std::max)(maxval, inPtr[traversal_info[in.shape().size()]]);
      traverse_next(in, shape_multipliers, traversal_info,
                    counter_offset + c + 1);
    }
    outPtr[r] = maxval;
  }
}

// =========Matrix operations ================================================
/*
template <>
void SoftMax<float, lang::Cpp>(const Tensor &in, Tensor *out, Context* ctx) {
  CHECK_LE(in.nDim(), 2u) << "Axis is required for SoftMax on multi dimemsional
tensor";
  out->CopyData(in);
  size_t nrow = 1, ncol = in.Size(), size = ncol;
  if (in.nDim() == 2u) {
    nrow = in.shape(0);
    ncol = size / nrow;
    out->Reshape(Shape{nrow, ncol});
  }
  Tensor tmp = RowMax(*out);
  SubColumn(tmp, out);
  Exp(*out, out);

  SumColumns(*out, &tmp);
  DivColumn(tmp, out);
  out->Reshape(in.shape());
}

template <>
void AddCol<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                              const Tensor& A, const Tensor& v, Tensor* out,
                              Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *APtr = static_cast<const float *>(A.data());
  const float *vPtr = static_cast<const float *>(v.data());
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    for (size_t c = 0; c < ncol; c++) {
      outPtr[offset + c] = APtr[offset + c] + vPtr[r];
    }
  }
}

template <>
void AddRow<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                              const Tensor& A, const Tensor& v, Tensor* out,
                              Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *APtr = static_cast<const float *>(A.data());
  const float *vPtr = static_cast<const float *>(v.data());
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    for (size_t c = 0; c < ncol; c++) {
      outPtr[offset + c] = APtr[offset + c] + vPtr[c];
    }
  }
}
template <>
void Outer<float, lang::Cpp>(const size_t m, const size_t n, const Tensor& in1,
                             const Tensor& in2, Tensor* out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *in1Ptr = static_cast<const float *>(in1.data());
  const float *in2Ptr = static_cast<const float *>(in2.data());
  for (size_t r = 0; r < m; r++) {
    size_t offset = r * n;
    for (size_t c = 0; c < n; c++) {
      outPtr[offset + c] = in1Ptr[r] * in2Ptr[c];
    }
  }
}
template <>
void Softmax<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                               const Tensor& in, Tensor* out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in.data());
  float *bPtr = new float[ncol];
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    float denom = 0.f;
    for (size_t c = 0; c < ncol; c++) {
      bPtr[c] = exp(inPtr[offset + c]);
      denom += bPtr[c];
    }
    for (size_t c = 0; c < ncol; c++) {
      size_t idx = offset + c;
      outPtr[idx] = bPtr[c] / denom;
    }
  }
  delete bPtr;
}

template <>
void SumColumns<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                                  const Tensor& in, Tensor* out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in.data());
  for (size_t c = 0; c < ncol; c++) {
    outPtr[c] = 0.f;
  }
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    for (size_t c = 0; c < ncol; c++) {
      outPtr[c] += inPtr[offset + c];
    }
  }
}

template <>
void SumRows<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                               const Tensor& in, Tensor* out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in.data());
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    outPtr[r] = 0.f;
    for (size_t c = 0; c < ncol; c++) {
      outPtr[r] += inPtr[offset + c];
    }
  }
}
*/
}  // namespace singa

#endif  // SINGA_CORE_TENSOR_TENSOR_MATH_CPP_H_

ext: h

max_stars_repo_name: fukien/incubator-singa

max_issues_repo_name: fukien/incubator-singa

max_forks_repo_name: fukien/incubator-singa


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <mpi.h>

/* Pull in the CBLAS header. */
#ifdef _MACOSX
#include <Accelerate/Accelerate.h>
#else
#ifdef _ATLAS
#include <cblas.h>
#else
#include <gsl_cblas.h>
#endif /* _ATLAS */
#endif /* _MACOSX */

/* These headers are provided by ScaleME. */
#include "ScaleME.h"

#include "precision.h"

#include "mlfma.h"
#include "direct.h"
#include "itsolver.h"
#include "util.h"

int matvec (cplx *out, cplx *in, cplx *cur, int id) {
	long i, nelt = (long)fmaconf.numbases * (long)fmaconf.bspboxvol;

	/* Compute the contrast pressure. */
#pragma omp parallel for default(shared) private(i)
	for (i = 0; i < nelt; ++i) cur[i] = in[i] * fmaconf.contrast[i];

	/* Reset the direct-interaction buffer and compute
	 * the matrix-vector product for the Green's matrix. */
	clrdircache();
	ScaleME_applyParFMA (cur, out);

	if (!id) return 0;

	/* Add in the identity portion. */
#pragma omp parallel for default(shared) private(i)
	for (i = 0; i < nelt; ++i) out[i] = fmaconf.cellvol * in[i] - out[i];

	return 0;
}

int gmres (cplx *rhs, cplx *sol, int guess,
		int mit, real tol, int quiet, augspace *aug) {
	long j, nelt = (long)fmaconf.numbases * (long)fmaconf.bspboxvol, lwork;
	int i, rank, one = 1, mred = mit;
	cplx *h, *v, *mvp, *beta, *vp, *hp, *s, cr,
		cone = 1., czero = 0., *azp, *zp;
	real rhn, err, *c;

	MPI_Comm_rank (MPI_COMM_WORLD, &rank);

	/* Allocate space for all required complex vectors. */
	lwork = (mit + 1) * (mit + nelt + 1) + nelt + mit;
	v = calloc (lwork, sizeof(cplx));	/* The Krylov subspace. */
	beta = v + nelt * (mit + 1);		/* The least-squares RHS. */
	mvp = beta + mit + 1;			/* Buffer for matrix-vector product. */
	h = mvp + nelt;				/* The upper Hessenberg matrix. */
	s = h + (mit + 1) * mit;		/* Givens rotation sines. */

	/* Allocate space for the Givens rotation cosines. */
	c = malloc (mit * sizeof(real));

	/* Compute the norm of the RHS for residual scaling. */
	rhn = parnorm(rhs, nelt);

	/* Compute the initial matrix-vector product for the input guess. */
	if (guess) matvec (v, sol, mvp, 1);

	/* Subtract from the RHS to form the residual. */
#pragma omp parallel for default(shared) private(j)
	for (j = 0; j < nelt; ++j) v[j] = rhs[j] - v[j];

	/* Zero the initial guess if one wasn't provided. */
	if (!guess) memset (sol, 0, nelt * sizeof(cplx));

	/* Find the norm of the initial residual. */
	err = parnorm(v, nelt);

	/* Construct the initial Arnoldi vector by normalizing the residual. */
#pragma omp parallel for default(shared) private(j)
	for (j = 0; j < nelt; ++j) v[j] /= err;

	/* Construct the vector beta for the minimization problem. */
	beta[0] = err;

	/* Report the RRE. */
	err /= rhn;
	if (!rank && !quiet) printf ("True residual: %g\n", err);

	/* The reduced number of iterations, if the space is augmented. */
	if (aug) mred = mit - aug->ntot;

	for (i = 0; i < mit && err > tol; ++i) {
		/* Point to the working space for this iteration. */
		vp = v + i * nelt;
		hp = h + i * (mit + 1);

		/* Compute the next expansion of the Krylov space. */
		if (!aug || i < mred) matvec (vp + nelt, vp, mvp, 1);
		else {
			/* Update with the next augmented vector. */
			azp = aug->az + nelt *
				((aug->nmax + aug->start + mred - i) % aug->nmax);
			/* Use the augmented space. */
			memcpy (vp + nelt, azp, nelt * sizeof(cplx));
		}

		/* Perform modified Gram-Schmidt to orthogonalize the basis. */
		/* This also builds the Hessenberg matrix column, including
		 * the 2-norm of the next basis vector. */
		cmgs (vp + nelt, hp, v, nelt, i + 1);

		/* Watch for breakdown. */
		if (cabs(hp[i + 1]) < REAL_EPSILON) {
			++i;
			break;
		}

		/* Apply previous Givens rotations to the Hessenberg column. */
		for (j = 0; j < i; ++j)
			ROT (&one, hp + j, &one, hp + j + 1, &one, c + j, s + j);

		/* Compute the Givens rotation for the current iteration. */
		LARTG (hp + i, hp + i + 1, c + i, s + i, &cr);
		/* Apply the current Givens rotation to the Hessenberg column. */
		hp[i] = cr;
		hp[i + 1] = 0;
		/* Perform the rotation on the vector beta. */
		ROT (&one, beta + i, &one, beta + i + 1, &one, c + i, s + i);

		/* Estimate the RRE for this iteration. */
		err = cabs(beta[i + 1]) / rhn;
		if (!rank && !quiet) printf ("GMRES(%d): %g\n", i, err);

		/* Flush the output buffers. */
		fflush (stdout);
		fflush (stderr);
	}

	/* If there were any GMRES iterations, update the solution. */
	if (i > 0 && aug) {
		/* Compute the optimum solution in the Krylov basis. */
		cplx *ys;
		ys = malloc(i * sizeof(cplx));
		memcpy (ys, beta, i * sizeof(cplx));
		TRSV (CblasColMajor, CblasUpper, CblasNoTrans,
				CblasNonUnit, i, h, mit + 1, ys, 1);

		/* Compute the next Krylov vector in the augmented space. */
		aug->start = (aug->start + 1) % aug->nmax;
		if (aug->ntot < aug->nmax) ++(aug->ntot);
		azp = aug->az + nelt * aug->start;

		beta[i] = 0;
		for (j = i - 1; j >= 0; --j) {
			/* Invert the Givens rotations. */
			s[j] = -s[j];
			ROT (&one, beta + j, &one, beta + j + 1, &one, c + j, s + j);
		}
		GEMV (CblasColMajor, CblasNoTrans, nelt, i + 1,
				&cone, v, nelt, beta, 1, &czero, azp, 1);

		/* Overwrite the Krylov subspace with the augmented subspace.
		 * The start pointer has already been incremented! */
		for (j = mred; j < i; ++j) {
			zp = aug->z + nelt *
				((aug->nmax + aug->start + mred - j - 1) % aug->nmax);
			memcpy(v + j * nelt, zp, nelt * sizeof(cplx));
		}

		/* Compute the solution update. */
		zp = aug->z + nelt * aug->start;
		GEMV (CblasColMajor, CblasNoTrans, nelt, i,
				&cone, v, nelt, ys, 1, &czero, zp, 1);
		for (j = 0; j < nelt; ++j) sol[j] += zp[j];

		free(ys);
	} else if (i > 0) {
		/* Compute the minimizer of the least-squares problem. */
		TRSV (CblasColMajor, CblasUpper, CblasNoTrans,
				CblasNonUnit, i, h, mit + 1, beta, 1);

		/* Compute the solution update in place. */
		GEMV (CblasColMajor, CblasNoTrans, nelt, i,
				&cone, v, nelt, beta, 1, &cone, sol, 1);
	}

	free (v);
	free (c);

	return i;
}

int bicgstab (cplx *rhs, cplx *sol,
		int guess, int mit, real tol, int quiet) {
	long j, nelt = (long)fmaconf.numbases * (long)fmaconf.bspboxvol;
	int i, rank;
	cplx *r, *rhat, *v, *p, *mvp, *t;
	cplx rho, alpha, omega, beta;
	real err, rhn;

	MPI_Comm_rank (MPI_COMM_WORLD, &rank);

	rho = alpha = omega = 1.;

	/* Allocate and zero the work arrays. */
	r = calloc (6L * nelt, sizeof(cplx));
	rhat = r + nelt;
	v = rhat + nelt;
	p = v + nelt;
	t = p + nelt;
	mvp = t + nelt;

	/* Compute the norm of the right-hand side for residual scaling. */
	rhn = parnorm(rhs, nelt);

	/* Compute the inital matrix-vector product for the input guess. */
	if (guess) matvec (r, sol, mvp, 1);

	/* Subtract from the RHS to form the residual. */
#pragma omp parallel for default(shared) private(j)
	for (j = 0; j < nelt; ++j) r[j] = rhs[j] - r[j];

	if (!guess) memset (sol, 0, nelt * sizeof(cplx));

	/* Copy the initial residual as the test vector. */
	memcpy (rhat, r, nelt * sizeof(cplx));

	/* Find the norm of the initial residual. */
	err = parnorm(r, nelt) / rhn;
	if (!rank && !quiet) printf ("True residual: %g\n", err);

	/* Run iterations until convergence or the maximum is reached. */
	for (i = 0; i < mit && err > tol; ++i) {
		/* Pre-compute portion of beta from previous iteration. */
		beta = alpha / (rho * omega);
		/* Compute rho for this iteration. */
		rho = pardot (rhat, r, nelt);
		/* Include the missing factor in beta. */
		beta *= rho;

		/* Update the search vector. */
#pragma omp parallel for default(shared) private(j)
		for (j = 0; j < nelt; ++j)
			p[j] = r[j] + beta * (p[j] - omega * v[j]);

		/* Compute the first search step, v = A * p. */
		matvec (v, p, mvp, 1);

		/* Compute the next alpha. */
		alpha = rho / pardot (rhat, v, nelt);

#pragma omp parallel for default(shared) private(j)
		for (j = 0; j < nelt; ++j) {
			/* Update the solution vector. */
			sol[j] += alpha * p[j];
			/* Update the residual vector. */
			r[j] -= alpha * v[j];
		}

		/* Compute the scaled residual norm and stop if convergence
		 * has been achieved. */
		err = parnorm(r, nelt) / rhn;
		if (!rank && !quiet) printf ("BiCG-STAB(%0.1f): %g\n", 0.5 + i, err);

		/* Flush the output buffers. */
		fflush (stdout);
		fflush (stderr);

		if (err < tol) break;

		/* Compute the next search step, t = A * r. */
		matvec (t, r, mvp, 1);

		/* Compute the update direction. */
		omega = pardot (t, r, nelt) / pardot (t, t, nelt);

		/* Update both the residual and the solution guess. */
#pragma omp parallel for default(shared) private(j)
		for (j = 0; j < nelt; ++j) {
			/* Update the solution vector. */
			sol[j] += omega * r[j];
			/* Update the residual vector. */
			r[j] -= omega * t[j];
		}

		/* Compute the scaled residual norm. */
		err = parnorm(r, nelt) / rhn;
		if (!rank && !quiet) printf ("BiCG-STAB(%d): %g\n", i + 1, err);

		fflush (stdout);
		fflush (stderr);
	}

	free (r);
	return i;
}

max_stars_repo_name: ahesford/afma

max_issues_repo_name: ahesford/afma

max_forks_repo_name: ahesford/afma


text:
#ifndef UTILS_H
#define UTILS_H

#ifdef __cplusplus
extern "C" {
#endif

#include <stdio.h> // for printf
#include <stdlib.h> // for posix_memalign
#include <string.h> // for memset
#include <math.h> // for sqrt
#include <float.h> //for FLT_MAX
#include <time.h> // for clock_gettime
#include <cblas.h>

#define STACK_SIZE_TOTAL 0x8000000 // 128Mbytes

typedef struct
    {
    int index;
    float score;
    }score_index_struct;

typedef struct
    {
    float x;
    float y;
    float w;
    float h;
    }box_struct;

typedef struct
    {
    void *stack_starting_address;
    char *stack_current_address;
    unsigned int stack_current_alloc_size;
    }stack_struct;

void init_stack
    (
    void
    );

void free_stack
    (
    void
    );

void *alloc_from_stack
    (
    unsigned int len
    );

void partial_free_from_stack
    (
    unsigned int len
    );

unsigned int get_stack_current_alloc_size
    (
    void
    );

void reset_stack_ptr_to_assigned_position
    (
    unsigned int assigned_size
    );

double what_time_is_it_now
    (
    void
    );

#ifdef __cplusplus
}
#endif

#endif // UTILS_H

max_stars_repo_name: lincolnhard/mobilenetssd-C-scratch

max_issues_repo_name: lincolnhard/mobilenetssd-C-scratch

max_forks_repo_name: lincolnhard/mobilenetssd-C-scratch


text:
#pragma once

#include <gsl/gsl>

#include "halley/file_formats/config_file.h"

namespace YAML {
	class Node;
	class Emitter;
}

namespace Halley {
	class Path;

	class YAMLConvert {
	public:
		class EmitOptions {
		public:
			Vector<String> mapKeyOrder;
		};

		static void parseConfig(ConfigFile& config, gsl::span<const gsl::byte> data);
		static ConfigNode parseYAMLNode(const YAML::Node& node);
		static ConfigFile parseConfig(gsl::span<const gsl::byte> data);
		static ConfigFile parseConfig(const Bytes& data);
		static ConfigNode parseConfig(const String& str);

		static String generateYAML(const ConfigFile& config, const EmitOptions& options);
		static String generateYAML(const ConfigNode& node, const EmitOptions& options);

	private:
		static void emitNode(const ConfigNode& node, YAML::Emitter& emitter, const EmitOptions& options);
		static void emitSequence(const ConfigNode& node, YAML::Emitter& emitter, const EmitOptions& options);
		static void emitMap(const ConfigNode& node, YAML::Emitter& emitter, const EmitOptions& options);
		static bool isCompactSequence(const ConfigNode& node, int depth);
	};

}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


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

#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <qdm.h>

gsl_vector *
qdm_knots_vector(size_t spline_df, gsl_vector *interior_knots)
{
  size_t size = spline_df * 2 + interior_knots->size;
  gsl_vector *knots = gsl_vector_alloc(size);

  // Fill front with zeros.
  for (size_t i = 0; i < spline_df; i++) {
    gsl_vector_set(knots, i, 0);
  }

  // Fill middle with interior knots.
  gsl_vector_view view = gsl_vector_subvector(knots, spline_df, interior_knots->size);
  gsl_blas_dcopy(interior_knots, &view.vector);

  // Fill end with ones.
  for (size_t i = spline_df + interior_knots->size; i < knots->size; i++) {
    gsl_vector_set(knots, i, 1);
  }

  return knots;
}

int
qdm_knots_rss(
  double *rss,

  gsl_vector *sorted_data,
  gsl_vector *middle,
  gsl_vector *interior_knots,

  size_t spline_df
)
{
  int status = 0;

  gsl_vector *m_knots = qdm_knots_vector(spline_df, interior_knots);

  size_t m = spline_df + interior_knots->size;
  gsl_matrix *ix = gsl_matrix_alloc(middle->size, m + 1);

  gsl_vector *emperical_quantiles = qdm_vector_quantile(sorted_data, middle);
  gsl_vector *theta = gsl_vector_alloc(ix->size2);
  gsl_vector *spline_quantiles = gsl_vector_alloc(emperical_quantiles->size);

  /* Build the true quantile function. */
  qdm_ispline_matrix(ix, middle, spline_df, m_knots);

  /* Check if the matrix is full rank. */
  double det = 0;
  status = qdm_matrix_det_tmm(ix, &det);
  if (status != 0) {
    goto cleanup;
  }

  if (det <= 0) {
    goto cleanup;
  }

  status = qdm_theta_optimize(theta, emperical_quantiles, ix);
  if (status != 0) {
    goto cleanup;
  }

  status = gsl_blas_dgemv(CblasNoTrans , 1.0, ix, theta, 0.0, spline_quantiles);
  if (status != 0) {
    goto cleanup;
  }

  *rss = qdm_vector_rss(emperical_quantiles, spline_quantiles);

cleanup:
  gsl_vector_free(spline_quantiles);
  gsl_vector_free(theta);
  gsl_vector_free(emperical_quantiles);
  gsl_matrix_free(ix);
  gsl_vector_free(m_knots);

  return status;
}

int
qdm_knots_optimize(
  gsl_vector *result,

  gsl_rng *rng,

  gsl_vector *sorted_data,
  gsl_vector *middle,
  gsl_vector *possible_knots,

  size_t iterate_n,
  size_t spline_df
)
{
  int status = 0;
  double rss = INFINITY;
  double rss_found = 0;

  gsl_vector *interior_knots = gsl_vector_alloc(result->size);

  for (size_t i = 0; i < iterate_n; i++) {
    status = gsl_ran_choose(
        rng,
        interior_knots->data,
        interior_knots->size,
        possible_knots->data,
        possible_knots->size,
        sizeof(double)
    );
    if (status != 0) {
      goto cleanup;
    }

    status = qdm_knots_rss(
        &rss_found,

        sorted_data,
        middle,
        interior_knots,

        spline_df
    );
    if (status != 0) {
      goto cleanup;
    }

    if (rss_found < rss) {
      rss = rss_found;

      status = gsl_vector_memcpy(result, interior_knots);
      if (status != 0) {
        goto cleanup;
      }
    }
  }

cleanup:
  gsl_vector_free(interior_knots);

  return status;
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
/**
 *
 * @file core_cgelqt.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:45 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_cgelqt - computes a LQ factorization of a complex M-by-N tile A: A = L * Q.
 *
 *  The tile Q is represented as a product of elementary reflectors
 *
 *    Q = H(k)' . . . H(2)' H(1)', where k = min(M,N).
 *
 *  Each H(i) has the form
 *
 *    H(i) = I - tau * v * v'
 *
 *  where tau is a complex scalar, and v is a complex vector with
 *  v(1:i-1) = 0 and v(i) = 1; conjfg(v(i+1:n)) is stored on exit in
 *  A(i,i+1:n), and tau in TAU(i).
 *
 *******************************************************************************
 *
 * @param[in] M
 *          The number of rows of the tile A.  M >= 0.
 *
 * @param[in] N
 *         The number of columns of the tile A.  N >= 0.
 *
 * @param[in] IB
 *         The inner-blocking size.  IB >= 0.
 *
 * @param[in,out] A
 *         On entry, the M-by-N tile A.
 *         On exit, the elements on and below the diagonal of the array
 *         contain the M-by-min(M,N) lower trapezoidal tile L (L is
 *         lower triangular if M <= N); the elements above the diagonal,
 *         with the array TAU, represent the unitary tile Q as a
 *         product of elementary reflectors (see Further Details).
 *
 * @param[in] LDA
 *         The leading dimension of the array A.  LDA >= max(1,M).
 *
 * @param[out] T
 *         The IB-by-N triangular factor T of the block reflector.
 *         T is upper triangular by block (economic storage);
 *         The rest of the array is not referenced.
 *
 * @param[in] LDT
 *         The leading dimension of the array T. LDT >= IB.
 *
 * @param[out] TAU
 *         The scalar factors of the elementary reflectors (see Further
 *         Details).
 *
 * @param[out] WORK
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgelqt = PCORE_cgelqt
#define CORE_cgelqt PCORE_cgelqt
#endif
int CORE_cgelqt(int M, int N, int IB,
                PLASMA_Complex32_t *A, int LDA,
                PLASMA_Complex32_t *T, int LDT,
                PLASMA_Complex32_t *TAU,
                PLASMA_Complex32_t *WORK)
{
    int i, k, sb;

    /* Check input arguments */
    if (M < 0) {
        coreblas_error(1, "Illegal value of M");
        return -1;
    }
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if ((IB < 0) || ( (IB == 0) && ((M > 0) && (N > 0)) )) {
        coreblas_error(3, "Illegal value of IB");
        return -3;
    }
    if ((LDA < max(1,M)) && (M > 0)) {
        coreblas_error(5, "Illegal value of LDA");
        return -5;
    }
    if ((LDT < max(1,IB)) && (IB > 0)) {
        coreblas_error(7, "Illegal value of LDT");
        return -7;
    }

    /* Quick return */
    if ((M == 0) || (N == 0) || (IB == 0))
        return PLASMA_SUCCESS;

    k = min(M, N);

    for(i = 0; i < k; i += IB) {
        sb = min(IB, k-i);

        LAPACKE_cgelq2_work(LAPACK_COL_MAJOR, sb, N-i,
                            &A[LDA*i+i], LDA, &TAU[i], WORK);

        LAPACKE_clarft_work(LAPACK_COL_MAJOR,
            lapack_const(PlasmaForward),
            lapack_const(PlasmaRowwise),
            N-i, sb,
            &A[LDA*i+i], LDA, &TAU[i],
            &T[LDT*i], LDT);

        if (M > i+sb) {
            LAPACKE_clarfb_work(
                LAPACK_COL_MAJOR,
                lapack_const(PlasmaRight),
                lapack_const(PlasmaNoTrans),
                lapack_const(PlasmaForward),
                lapack_const(PlasmaRowwise),
                M-i-sb, N-i, sb,
                &A[LDA*i+i],      LDA,
                &T[LDT*i],        LDT,
                &A[LDA*i+(i+sb)], LDA,
                WORK, M-i-sb);
        }
    }
    return PLASMA_SUCCESS;
}


text:
#pragma once
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

double integrator_qag(gsl_function *, double, double, double, double, size_t, double *,int *);
double integrator_cquad(gsl_function *,double,double,double,double,size_t,double *,int *);
double integrator_qags(gsl_function *,double,double,double,double,size_t,double *,int *);


text:
/* Copyright (c) 2017 Bradley Worley <geekysuavo@gmail.com>
 * Released under the MIT License
 */

/* ensure once-only inclusion. */
#ifndef __MATTE_BLAS_H__
#define __MATTE_BLAS_H__

/* include the vector and matrix headers. */
#include <matte/vector.h>
#include <matte/matrix.h>

/* include the complex vector and matrix headers. */
#include <matte/complex-vector.h>
#include <matte/complex-matrix.h>

/* include the atlas blas header. */
#include <cblas.h>

/* function declarations, level 1, double (blas.c): */

int matte_daxpy (double alpha, Vector x, Vector y);

int matte_dscal (double alpha, Vector x);

int matte_dcopy (Vector x, Vector y);

int matte_dswap (Vector x, Vector y);

int matte_ddot (Vector x, Vector y, double *out);

int matte_dnrm2 (Vector x, double *out);

/* function declarations, level 1, complex double (blas.c): */

int matte_zaxpy (complex double alpha, ComplexVector x, ComplexVector y);

int matte_zscal (complex double alpha, ComplexVector x);

int matte_zcopy (ComplexVector x, ComplexVector y);

int matte_zswap (ComplexVector x, ComplexVector y);

int matte_zdotu (ComplexVector x, ComplexVector y, complex double *out);

int matte_zdotc (ComplexVector x, ComplexVector y, complex double *out);

int matte_znrm2 (ComplexVector x, double *out);

/* function declarations, level 2, double (blas.c): */

int matte_dgemv (MatteTranspose trans, double alpha, Matrix A, Vector x,
                 double beta, Vector y);

int matte_dsymv (MatteTriangle uplo, double alpha, Matrix A, Vector x,
                 double beta, Vector y);

int matte_dtrmv (MatteTriangle uplo, MatteTranspose trans,
                 MatteDiagonal diag, Matrix A, Vector x);

int matte_dtrsv (MatteTriangle uplo, MatteTranspose trans,
                 MatteDiagonal diag, Matrix A, Vector x);

int matte_dger (double alpha, Vector x, Vector y, Matrix A);

int matte_dsyr (MatteTriangle uplo, double alpha, Vector x, Matrix A);

int matte_dsyr2 (MatteTriangle uplo, double alpha, Vector x, Vector y,
                 Matrix A);

/* function declarations, level 2, complex double (blas.c): */

int matte_zgemv (MatteTranspose trans, complex double alpha, ComplexMatrix A,
                 ComplexVector x, complex double beta, ComplexVector y);

int matte_zhemv (MatteTriangle uplo, complex double alpha, ComplexMatrix A,
                 ComplexVector x, complex double beta, ComplexVector y);

int matte_ztrmv (MatteTriangle uplo, MatteTranspose trans,
                 MatteDiagonal diag, ComplexMatrix A,
                 ComplexVector x);

int matte_ztrsv (MatteTriangle uplo, MatteTranspose trans,
                 MatteDiagonal diag, ComplexMatrix A,
                 ComplexVector x);

int matte_zher (MatteTriangle uplo, double alpha,
                ComplexVector x, ComplexMatrix A);

int matte_zher2 (MatteTriangle uplo, complex double alpha,
                 ComplexVector x, ComplexVector y,
                 ComplexMatrix A);

/* function declarations, level 3, double (blas.c): */

int matte_dgemm (MatteTranspose transA, MatteTranspose transB,
                 double alpha, Matrix A, Matrix B,
                 double beta, Matrix C);

/* FIXME */

/* function declarations, level 3, complex double (blas.c): */

int matte_zgemm (MatteTranspose transA, MatteTranspose transB,
                 complex double alpha, ComplexMatrix A, ComplexMatrix B,
                 complex double beta, ComplexMatrix C);

/* FIXME */

#endif /* !__MATTE_BLAS_H__ */

max_stars_repo_name: geekysuavo/matte

max_issues_repo_name: geekysuavo/matte

max_forks_repo_name: geekysuavo/matte


text:
#include <cblas.h>

#include "tasks.h"

void solve_task_seq(void *ptr)
{
    struct solve_task_arg *arg = (struct solve_task_arg*) ptr;

    int n     = arg->n;
    int nrhs  = arg->nrhs;
    double *L = arg->L;
    int ldL   = arg->ldL;
    double *X = arg->X;
    int ldX   = arg->ldX;

    // Solve L * X = B, where B and X share the same memory.
    cblas_dtrsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit,
                n, nrhs,
                1.0, L, ldL,
                     X, ldX);
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef SIM_PHYSICS_ROTOR_DATABASE_3D_H_
#define SIM_PHYSICS_ROTOR_DATABASE_3D_H_

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

#include <string>

#include "common/macros.h"

class RotorDatabase3d {
 public:
  explicit RotorDatabase3d(const std::string &filename);
  virtual ~RotorDatabase3d();

  // The rotor force prime frame is attached to an isolated rotor, and is
  // defined as follows:
  // X_prime: The rotor shaft axis, positive in the thrusting direction.
  // Y_prime: The in-plane axis in the Z_prime cross X_prime direction.
  // Z_prime: The in-plane axis aligned with the apparent edgewise wind
  //          direction. Positive is into the wind.

  // Z_prime force [N] is a rotor in-plane force [N].
  double CalcForceZPrime(double angular_rate, double axial_vel,
                         double edgewise_vel, double air_density) const;

  // Y_prime force [N] is the rotor in-plane side force [N].
  // The sign is flipped if rotor direction does not match direction the
  // the database was created assuming.
  double CalcForceYPrime(double angular_rate, double axial_vel,
                         double edgewise_vel, double air_density,
                         int dir) const;

  // Roll moment [N-m] is the right-handed moment about the Z_prime force axis.
  double CalcMomZPrime(double angular_rate, double axial_vel,
                       double edgewise_vel, double air_density, int dir) const;

  // Pitch moment [N-m] is the right-handed moment about the Y_prime axis:
  // positive is pitch up/back.
  double CalcMomYPrime(double angular_rate, double axial_vel,
                       double edgewise_vel, double air_density) const;

  // Thrust [N] is positive for propulsive/motoring, and is equivalent
  // to the rotor X or X_prime axis force.
  double CalcThrust(double angular_rate, double axial_vel, double edgewise_vel,
                    double air_density) const;

  // Torque [N-m] is aerodynamic moment about thrust axis (opposes rotation
  // direction for positive thrusts).
  double CalcTorque(double angular_rate, double axial_vel, double edgewise_vel,
                    double air_density, int dir) const;

  double LookupZForceCoeff(double angular_rate, double axial_vel,
                           double edgewise_vel) const;
  double LookupYForceCoeff(double angular_rate, double axial_vel,
                           double edgewise_vel) const;
  double LookupZMomCoeff(double angular_rate, double axial_vel,
                         double edgewise_vel) const;
  double LookupYMomCoeff(double angular_rate, double axial_vel,
                         double edgewise_vel) const;
  double LookupThrustCoeff(double angular_rate, double axial_vel,
                           double edgewise_vel) const;
  double LookupTorqueCoeff(double angular_rate, double axial_vel,
                           double edgewise_vel) const;
  double CalcThrustScaleFactor(double angular_rate, double axial_vel,
                               double edgewise_vel) const;
  double CalcTorqueScaleFactor(double angular_rate, double axial_vel,
                               double edgewise_vel) const;

 private:
  double radius_;
  double radius_pow4_;
  double radius_pow5_;
  int32_t spin_;

  gsl_vector *angular_rates_;
  gsl_vector *axial_vels_;
  gsl_vector *edgewise_vels_;

  gsl_vector *z_prime_force_coeffs_;
  gsl_vector *y_prime_force_coeffs_;
  gsl_vector *z_prime_mom_coeffs_;
  gsl_vector *y_prime_mom_coeffs_;
  gsl_vector *thrust_coeffs_;
  gsl_vector *torque_coeffs_;

  DISALLOW_COPY_AND_ASSIGN(RotorDatabase3d);
};

#endif  // SIM_PHYSICS_ROTOR_DATABASE_3D_H_


text:
// CONSTANTIN MIHAI - 336CA
// TEMA 2 ASC

#include <cblas.h>

#include "utils.h"

double* copy(int N, double *A) {

	int i, j;

	double *B = malloc(N * N * sizeof(double));
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			B[i * N + j] = A[i * N + j];
		}
	}

	return B;
}

double* add_matrix(int N, double *A, double *B) {

	int i, j;

	double *C = calloc(N * N, sizeof(double));
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			C[i * N + j] = A[i * N + j] + B[i * N + j];
		}
	}
	return C;
}

double* my_solver(int N, double *A, double *B) {

	/*
		A - matrice superior triunghiulara
		B - matrice normala
		Vrem C = B * A^t + A^2 * B
	*/

	/*
		CblasRowMajor - matrice stocata pe o singura linie
		CblasRight - facem B * A^t
		CblasUpper - matrice superior triunghiulara
		CblasTrans - folosim A^t
		CblasNoTrans - folosim A, nu A^t
		CblasNonUnit - A nu e matrice unitara
	*/

	// B2 <- B2 * A^t
	double *B2 = copy(N, B);
	cblas_dtrmm(CblasRowMajor, CblasRight, CblasUpper, CblasTrans, CblasNonUnit, N, N, 1.0, A, N, B2, N);

	// A2 <- A2 * A
	double *A2 = copy(N, A);
	cblas_dtrmm(CblasRowMajor, CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, N, N, 1.0, A, N, A2, N);

	// B <- A2 * B
	cblas_dtrmm(CblasRowMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, N, N, 1.0, A2, N, B, N);

	// result <- B2 + B
	double *result = add_matrix(N, B2, B);

	return result;
}

max_stars_repo_name: mihai-constantin/ACS

max_issues_repo_name: mihai-constantin/ACS

max_forks_repo_name: mihai-constantin/ACS