#include <stddef.h></stddef.h>

#include <stdlib.h></stdlib.h>

#include <memory></memory>

#include <cstring></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);</double>

	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>

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

	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);

}