#include <bench/refblas.h>
void
dgemv_ref(int m, int n, double alpha,
const double *A, int incRowA, int incColA,
const double *x, int incX,
double beta,
double *y, int incY)
{
int i, j;
for (i=0; i<m; ++i) {
if (beta==0) {
y[i*incY] = 0;
} else {
y[i*incY] *= beta;
}
for (j=0; j<n; ++j) {
y[i*incY] += alpha*A[i*incRowA+j*incColA]*x[j*incX];
}
}
}
void
dtrlmv_ref(int n, int unitDiag,
const double *A, int incRowA, int incColA,
double *x, int incX)
{
int i, j;
for (j=n-1; j>=0; --j) {
for (i=j+1; i<n; ++i) {
x[i*incX] += A[i*incRowA+j*incColA]*x[j*incX];
}
if (!unitDiag) {
x[j*incX] *= A[j*(incRowA+incColA)];
}
}
}
void
dtrlsv_ref(int n, int unitDiag,
const double *A, int incRowA, int incColA,
double *x, int incX)
{
int i, j;
for (i=0; i<n; ++i) {
for (j=0; j<i; ++j) {
x[i*incX] -= A[i*incRowA+j*incColA]*x[j*incX];
}
x[i*incX] /= (unitDiag) ? 1.0 : A[i*(incRowA+incColA)];
}
}
void
dger_ref(int m, int n, double alpha,
const double *x, int incX,
const double *y, int incY,
double *A, int incRowA, int incColA)
{
int i, j;
for (j=0; j<n; ++j) {
for (i=0; i<m; ++i) {
A[i*incRowA+j*incColA] += alpha*x[i*incX]*y[j*incY];
}
}
}
void
dgemm_ref(int m, int n, int k,
double alpha,
const double *A, int incRowA, int incColA,
const double *B, int incRowB, int incColB,
double beta,
double *C, int incRowC, int incColC)
{
int i, j, l;
if (beta!=1) {
for (i=0; i<m; ++i) {
for (j=0; j<n; ++j) {
if (beta!=0) {
C[i*incRowC+j*incColC] *= beta;
} else {
C[i*incRowC+j*incColC] = 0;
}
}
}
}
if (alpha!=0) {
for (i=0; i<m; ++i) {
for (j=0; j<n; ++j) {
for (l=0; l<k; ++l) {
C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
*B[l*incRowB+j*incColB];
}
}
}
}
}
void
dtrlsm_ref(int m, int n, double alpha, int unitDiag,
const double *A, int incRowA, int incColA,
double *B, int incRowB, int incColB)
{
int i, j;
if (alpha!=1) {
for (j=0; j<n; ++j) {
for (i=0; i<m; ++i) {
B[i*incRowB+j*incColB] *= alpha;
}
}
}
for (j=0; j<n; ++j) {
dtrlsv_ref(m, unitDiag, A, incRowA, incColA, &B[j*incColB], incRowB);
}
}