#include <flens/flens.cxx>
#ifndef MIN_M
#define MIN_M 50
#endif
#ifndef MIN_N
#define MIN_N 50
#endif
#ifndef MIN_K
#define MIN_K 50
#endif
#ifndef MAX_M
#define MAX_M 2000
#endif
#ifndef MAX_N
#define MAX_N 2000
#endif
#ifndef MAX_K
#define MAX_K 2000
#endif
#ifndef INC_M
#define INC_M 50
#endif
#ifndef INC_N
#define INC_N 50
#endif
#ifndef INC_K
#define INC_K 50
#endif
#ifndef ALPHA
#define ALPHA 1
#endif
#ifndef BETA
#define BETA 1
#endif
#ifndef TYPE
#define TYPE double
#endif
#ifndef TYPE_A
#define TYPE_A TYPE
#endif
#ifndef TYPE_B
#define TYPE_B TYPE
#endif
#ifndef TYPE_C
#define TYPE_C TYPE
#endif
#ifndef TYPE_X
#define TYPE_X TYPE
#endif
#ifndef TYPE_Y
#define TYPE_Y TYPE
#endif
#ifndef TYPE_ALPHA
#define TYPE_ALPHA TYPE
#endif
#ifndef TYPE_BETA
#define TYPE_BETA TYPE
#endif
template <typename T>
struct WallTime
{
void
tic()
{
t0 = std::chrono::high_resolution_clock::now();
}
T
toc()
{
using namespace std::chrono;
elapsed = high_resolution_clock::now() - t0;
return duration<T,seconds::period>(elapsed).count();
}
std::chrono::high_resolution_clock::time_point t0;
std::chrono::high_resolution_clock::duration elapsed;
};
template <typename T>
struct ComplexTrait
{
typedef T PT;
template <typename RD, typename DT>
static PT
rand(RD &rd, DT dt)
{
return dt(rd);
}
};
template <typename T>
struct ComplexTrait<std::complex<T> >
{
typedef T PT;
template <typename RD, typename DT>
static std::complex<T>
rand(RD &rd, DT dt)
{
return std::complex<T>(dt(rd), dt(rd));
}
};
// fill rectangular matrix with random values
template <typename MA>
void
fill(flens::GeMatrix<MA> &A)
{
typedef typename flens::GeMatrix<MA>::ElementType T;
typedef typename ComplexTrait<T>::PT PT;
typedef typename flens::GeMatrix<MA>::IndexType Index;
std::random_device random;
std::default_random_engine mt(random());
std::uniform_real_distribution<PT> uniform(-100,100);
for (Index i=1; i<=A.numRows(); ++i) {
for (Index j=1; j<=A.numCols(); ++j) {
A(i,j) = ComplexTrait<T>::rand(mt, uniform);
}
}
}
// fill dense vector with random values
template <typename VX>
void
fill(flens::DenseVector<VX> &x)
{
typedef typename flens::DenseVector<VX>::ElementType T;
typedef typename ComplexTrait<T>::PT PT;
typedef typename flens::DenseVector<VX>::IndexType Index;
std::random_device random;
std::default_random_engine mt(random());
std::uniform_real_distribution<PT> uniform(-100,100);
for (Index i=1; i<=x.length(); ++i) {
x(i) = ComplexTrait<T>::rand(mt, uniform);
}
}
// fill triangular matrix with random values and scale
// such that it gets diagonal dominant
template <typename MA>
void
fill(flens::TrMatrix<MA> &A)
{
typedef typename flens::TrMatrix<MA>::ElementType T;
typedef typename ComplexTrait<T>::PT PT;
typedef typename flens::TrMatrix<MA>::IndexType Index;
std::random_device random;
std::default_random_engine mt(random());
std::uniform_real_distribution<PT> uniform(-100,100);
PT thresh = 1.05;
if (A.upLo()==flens::Lower) {
for (Index i=1; i<=A.numRows(); ++i) {
PT rowSum = 0;
for (Index j=1; j<i; ++j) {
A(i,j) = ComplexTrait<T>::rand(mt, uniform);
rowSum += std::abs(A(i,j));
}
if (A.diag()!=flens::Unit) {
A(i,i) = ComplexTrait<T>::rand(mt, uniform);
if (std::abs(A(i,i))==PT(0)) {
A(i,i) = 1;
} else {
while (std::abs(A(i,i))<1) {
A(i,i) *= 2;
}
}
rowSum /= std::abs(A(i,i));
}
for (Index j=1; j<i; ++j) {
A(i,j) /= rowSum*thresh;
}
}
} else {
for (Index i=1; i<=A.numRows(); ++i) {
PT rowSum = 0;
for (Index j=i+1; j<A.numCols(); ++j) {
A(i,j) = ComplexTrait<T>::rand(mt, uniform);
rowSum += std::abs(A(i,j));
}
if (A.diag()!=flens::Unit) {
A(i,i) = ComplexTrait<T>::rand(mt, uniform);
if (std::abs(A(i,i))==PT(0)) {
A(i,i) = 1;
} else {
while (std::abs(A(i,i))<1) {
A(i,i) *= 2;
}
}
rowSum /= std::abs(A(i,i));
}
for (Index j=i+1; j<A.numCols(); ++j) {
A(i,j) /= rowSum*thresh;
}
}
}
}
template <typename MATRIX>
typename ComplexTrait<typename MATRIX::ElementType>::PT
asum(const MATRIX &A)
{
typedef typename MATRIX::ElementType T;
typedef typename ComplexTrait<T>::PT PT;
typedef typename MATRIX::IndexType Index;
PT asum = 0;
for (Index i=1; i<=A.numRows(); ++i) {
for (Index j=1; j<=A.numCols(); ++j) {
asum += std::abs(A(i,j));
}
}
return asum;
}
template <typename MA, typename MB>
typename ComplexTrait<typename MA::ElementType>::PT
asumDiff(const flens::GeMatrix<MA> &A, const flens::GeMatrix<MB> &B)
{
typedef typename flens::GeMatrix<MA>::ElementType T;
typedef typename ComplexTrait<T>::PT PT;
typedef typename flens::GeMatrix<MA>::IndexType Index;
PT asum = 0;
for (Index i=1; i<=A.numRows(); ++i) {
for (Index j=1; j<=A.numCols(); ++j) {
asum += std::abs(A(i,j) - B(i,j));
}
}
return asum;
}
template <typename VX, typename VY>
typename ComplexTrait<typename VX::ElementType>::PT
asumDiff(const flens::DenseVector<VX> &x, const flens::DenseVector<VY> &y)
{
typedef typename flens::DenseVector<VX>::ElementType T;
typedef typename ComplexTrait<T>::PT PT;
typedef typename flens::DenseVector<VX>::IndexType Index;
PT asum = 0;
for (Index i=1; i<=x.length(); ++i) {
asum += std::abs(x(i) - y(i));
}
return asum;
}
template <typename MA, typename MB, typename MC0, typename MC1>
double
estimateGemmResidual(const MA &A, const MB &B,
const MC0 &C0, const MC1 &C1)
{
typedef typename MC0::IndexType Index;
Index m= C1.numRows();
Index n= C1.numCols();
Index k= A.numCols();
double aNorm = asum(A);
double bNorm = asum(B);
double cNorm = asum(C1);
double diff = asumDiff(C0, C1);
// Using eps for double gives upper bound in case elements have lower
// precision.
double eps = std::numeric_limits<double>::epsilon();
double res = diff/(aNorm*bNorm*cNorm*eps*std::max(std::max(m,n),k));
return res;
}