#include <stddef.h>
#include <stdlib.h>
#include <memory>
#include <cstring>
 
//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);
}