#include <stddef.h>
#include <stdlib.h>
#include <memory>

#if defined(LINUX)
#include <cstring>
#endif

//blasint may either be common 4 bytes or extended 8 (long)...
//Replace this and REBUILD the CBLAS with extended int if needed.
typedef int blasint;

extern "C" {
#include "tlin/cblas.h"
}

#include "tlin/tlin_cblas_wrap.h"

//========================================================================================

/*
  LITTLE RESUME OF BLAS NOMENCLATURE

  Function names behave like:

    cblas_[Precision][Matrix type][Operation specifier]

  Precision can be:

    S - REAL                C - COMPLEX
    D - DOUBLE PRECISION    Z - COMPLEX*16  (may not be supported by all machines)

    In current wrapper implementation only D routines will be adopted.

  Matrix type can be:

    GE - GEneral     GB - General Band
    SY - SYmmetric   SB - Sym. Band     SP - Sym. Packed
    HE - HErmitian   HB - Herm. Band    HP - Herm. Packed
    TR - TRiangular  TB - Triang. Band  TP - Triang. Packed

    In current wrapper implementation only GE routines will be adopted.

  Operation specifier depends on the operation. For example, mm means matrix * matrix, mv matrix * vector, etc..
*/

//========================================================================================

void sum(int n, const double *x, double *&y)
{
	/*
    void cblas_daxpy(blasint n, double a, double *x, blasint incx, double *y, blasint incy);

    NOTE:
    
      Returns a * x + y.    Output is stored in y.

      incx and incy are the increments in array access - ie x[incx * i] and y[incy * j] values only are
      considered (=> we'll use 1). 
  */

	double *_x = const_cast<double *>(x);

	cblas_daxpy(n, 1.0, _x, 1, y, 1);
}

//---------------------------------------------------------------------------------------------------

void tlin::multiply(int rows, int cols, const double *A, const double *x, double *&y)
{
	/*
    void cblas_dgemv(enum CBLAS_ORDER order,  enum CBLAS_TRANSPOSE trans,  blasint m, blasint n,
		 double alpha, double  *a, blasint lda,  double  *x, blasint incx,  double beta,  double  *y, blasint incy);

     NOTE: 

      Returns alpha*A*x + beta*y.    Output is stored in y.
  */

	if (!y) {
		y = (double *)malloc(rows * sizeof(double));
		memset(y, 0, rows * sizeof(double));
	}

	double *_A = const_cast<double *>(A);
	double *_x = const_cast<double *>(x);

	cblas_dgemv(CblasColMajor, CblasNoTrans, rows, cols, 1.0, _A, rows, _x, 1, 1.0, y, 1);
}

//---------------------------------------------------------------------------------------------------

void tlin::multiply(const mat &A, const double *x, double *&y)
{
	multiply(A.rows(), A.cols(), A.values(), x, y);
}