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

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

#include <memory></memory>

#include <cstring></cstring>

#if defined(_MSC_VER) && (_MSC_VER >= 1900)

// dummy definition for linker

#include <conio.h></conio.h>

extern "C" {

void __imp__cprintf(char const *const _Format) { _cprintf(_Format); }

}

#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..

*/

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

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

}