7 #ifndef LIMBO_SOLVERS_NUMERICAL_H
8 #define LIMBO_SOLVERS_NUMERICAL_H
30 inline void mklaxpy(
unsigned int n, T a, T
const* x, T* y);
37 inline void mklaxpy<float>(
unsigned int n,
float a,
float const* x,
float* y)
39 cblas_saxpy(n, a, x, 1, y, 1);
47 inline void mklaxpy<double>(
unsigned int n,
double a,
double const* x,
double* y)
49 cblas_daxpy(n, a, x, 1, y, 1);
60 template <
typename T,
typename MatrixType>
61 inline void mklAxPlusy(T a, MatrixType
const& A, T
const* x, T b, T* y);
69 inline void mklAxPlusy<float, MatrixCSR<float, int, 1> >(
float a, MatrixCSR<float, int, 1>
const& A,
float const* x,
float b,
float* y)
71 mkl_scsrmv(
"N", &A.numRows, &A.numColumns, &a,
"G", A.vElement, A.vColumn, A.vRowBeginIndex, A.vRowBeginIndex+1, x, &b, y);
80 inline void mklAxPlusy<double, MatrixCSR<double, int, 1> >(
double a, MatrixCSR<double, int, 1>
const& A,
double const* x,
double b,
double* y)
82 mkl_dcsrmv(
"N", &A.numRows, &A.numColumns, &a,
"G", A.vElement, A.vColumn, A.vRowBeginIndex, A.vRowBeginIndex+1, x, &b, y);
93 template <
typename T,
typename MatrixType>
94 inline void mklATxPlusy(T a, MatrixType
const& A, T
const* x, T b, T* y);
101 inline void mklATxPlusy<float, MatrixCSR<float, int, 1> >(
float a, MatrixCSR<float, int, 1>
const& A,
float const* x,
float b,
float* y)
103 mkl_scsrmv(
"T", &A.numRows, &A.numColumns, &a,
"G", A.vElement, A.vColumn, A.vRowBeginIndex, A.vRowBeginIndex+1, x, &b, y);
111 inline void mklATxPlusy<double, MatrixCSR<double, int, 1> >(
double a, MatrixCSR<double, int, 1>
const& A,
double const* x,
double b,
double* y)
113 mkl_dcsrmv(
"T", &A.numRows, &A.numColumns, &a,
"G", A.vElement, A.vColumn, A.vRowBeginIndex, A.vRowBeginIndex+1, x, &b, y);
121 template <
typename T>
122 inline T mkldot(
unsigned int n, T
const* x, T
const* y);
129 inline float mkldot(
unsigned int n,
float const* x,
float const* y)
131 return cblas_sdot(n, x, 1, y, 1);
139 inline double mkldot(
unsigned int n,
double const* x,
double const* y)
141 return cblas_ddot(n, x, 1, y, 1);
149 template <
typename T>
150 inline void mklcopy(
unsigned int n, T
const* x, T* y);
156 inline void mklcopy(
unsigned int n,
float const* x,
float* y)
158 cblas_scopy(n, x, 1, y, 1);
165 inline void mklcopy(
unsigned int n,
double const* x,
double* y)
167 cblas_dcopy(n, x, 1, y, 1);
178 template <
typename T,
typename V>
179 inline void axpy(
unsigned int n, T a, V
const* x, T* y)
184 for (
unsigned int i = 0; i < n; ++i)
197 template <
typename T,
typename V,
typename MatrixType>
198 inline void AxPlusy(T a, MatrixType
const& A, V
const* x, T* y)
202 mklAxPlusy(a, A, x, (T)1, y);
204 for (
typename MatrixType::index_type i = 0; i < A.numRows; ++i)
206 for (
typename MatrixType::index_type k = A.vRowBeginIndex[i]; k < A.vRowBeginIndex[i+1]; ++k)
208 typename MatrixType::index_type kk = k-MatrixType::s_startingIndex;
209 #ifdef DEBUG_NUMERICAL
210 limboAssert(A.vColumn[kk]-MatrixType::s_startingIndex < A.numColumns && kk < A.numElements);
212 y[i] += a*A.vElement[kk]*x[A.vColumn[kk]-MatrixType::s_startingIndex];
225 template <
typename T,
typename V,
typename MatrixType>
226 inline void ATxPlusy(T a, MatrixType
const& A, V
const* x, T* y)
230 mklATxPlusy(a, A, x, (T)1, y);
232 for (
typename MatrixType::index_type i = 0; i < A.numRows; ++i)
234 for (
typename MatrixType::index_type k = A.vRowBeginIndex[i]; k < A.vRowBeginIndex[i+1]; ++k)
236 typename MatrixType::index_type kk = k-MatrixType::s_startingIndex;
237 #ifdef DEBUG_NUMERICAL
238 limboAssert(A.vColumn[kk]-MatrixType::s_startingIndex < A.numColumns && kk < A.numElements);
240 y[A.vColumn[kk]-MatrixType::s_startingIndex] += a*A.vElement[kk]*x[i];
251 template <
typename T>
252 inline T
dot(
unsigned int n, T
const* x, T
const* y)
256 result = mkldot(n, x, y);
258 for (
unsigned int i = 0; i < n; ++i)
268 template <
typename T>
269 inline void vcopy(
unsigned int n, T
const* x, T* y)
275 std::copy(x, x+n, y);
Basic utilities such as variables and linear expressions in solvers.
void ATxPlusy(T a, MatrixType const &A, V const *x, T *y)
void axpy(unsigned int n, T a, V const *x, T *y)
void vcopy(unsigned int n, T const *x, T *y)
copy vector
#define limboAssert(condition)
custom assertion without message
void AxPlusy(T a, MatrixType const &A, V const *x, T *y)
T dot(unsigned int n, T const *x, T const *y)
compute dot product