#include <catch2/catch.hpp>
#include <iostream>
#include <fenv.h>
#include <algorithm>
#include <limits>
#include <numeric>
#include "debug.h"
#include "mathfunc.h"
#include "fastblas.h"
#include "fastopenblas.h"
#include "param.h"
using namespace std;
TEST_CASE("calcPab", "[math]") {
REQUIRE(GetabIndex(1,1,3) == 0);
REQUIRE(GetabIndex(1,2,3) == 1);
REQUIRE(GetabIndex(1,3,3) == 2);
REQUIRE(GetabIndex(2,2,3) == 5);
REQUIRE(GetabIndex(2,3,3) == 6);
REQUIRE(GetabIndex(3,3,3) == 9);
REQUIRE(GetabIndex(2,1,3) == 1); // reverse
}
TEST_CASE( "Math functions", "[math]" ) {
debug_set_debug_mode(true);
debug_set_no_check_mode(false);
debug_set_strict_mode(true);
double data[] = { 2,-1, 0,
-1, 2,-1,
0,-1, 2};
gsl_matrix *m = gsl_matrix_alloc(3,3);
copy(data, data+9, m->data);
REQUIRE( isMatrixPositiveDefinite(m) );
REQUIRE( isMatrixSymmetric(m) );
// REQUIRE( checkMatrixEigen(m,0.001) );
double data1[] = {1.0,0,0,
0,3.0,0,
0,0,2.0};
copy(data1, data1+9, m->data);
REQUIRE( isMatrixPositiveDefinite(m) );
// REQUIRE( checkMatrixEigen(m) );
double data2[] = {1,1,1,
1,1,1,
1,1,0.5};
copy(data2, data2+9, m->data);
REQUIRE( !isMatrixPositiveDefinite(m));
// REQUIRE( !checkMatrixEigen(m) );
double data3[] = {1.0, 0, 0,
3.0,3.0, 0,
0, 0,2.0};
copy(data3, data3+9, m->data);
REQUIRE( !isMatrixPositiveDefinite(m) );
REQUIRE( !isMatrixSymmetric(m) );
// REQUIRE( checkMatrixEigen(m) );
// ---- NaN checks
vector<double> v = {1.0, 2.0};
// REQUIRE (!isnan(std::accumulate(v.begin(), v.end(), 0)));
vector<double> v2 = {1.0, 2.0, std::numeric_limits<double>::quiet_NaN()};
REQUIRE (isnan(v2[2]));
REQUIRE(has_nan(v2));
// test minus nan
vector<double> v3 = {1.0, 2.0, -std::numeric_limits<double>::quiet_NaN()};
REQUIRE (isnan(v3[2]));
REQUIRE(has_nan(v3));
}
TEST_CASE("cblas_dgemm", "[math]") {
double *A, *B, *C;
int m, n, k, i, j;
double alpha, beta;
printf ("\n This example computes real matrix C=alpha*A*B+beta*C using \n"
" Intel(R) MKL function dgemm, where A, B, and C are matrices and \n"
" alpha and beta are double precision scalars\n\n");
m = 2000, k = 200, n = 1000;
printf (" Initializing data for matrix multiplication C=A*B for matrix \n"
" A(%ix%i) and matrix B(%ix%i)\n\n", m, k, k, n);
alpha = 1.0; beta = 0.0;
printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n"
" performance \n\n");
A = (double *)malloc( m*k*sizeof( double ));
B = (double *)malloc( k*n*sizeof( double ));
C = (double *)malloc( m*n*sizeof( double ));
printf (" Intializing matrix data \n\n");
for (i = 0; i < (m*k); i++) {
A[i] = (double)(i+1);
}
for (i = 0; i < (k*n); i++) {
B[i] = (double)(-i-1);
}
for (i = 0; i < (m*n); i++) {
C[i] = 0.0;
}
printf (" Computing matrix product using Intel(R) MKL dgemm function via CBLAS interface \n\n");
assert(m==2000);
assert(k==200);
assert(n==1000);
//cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
// m, n, k, alpha, A, k, B, n, beta, C, n);
fast_cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
m, n, k, alpha, A, k, B, n, beta, C, n);
REQUIRE(trunc(C[0]) == -2666620100.0 );
REQUIRE(trunc(C[1]) == -2666640200.0 );
REQUIRE(trunc(C[2003]) == -10627000400.0 );
}
TEST_CASE("fast_dgemm", "[math]") {
double *A, *B, *C;
int m, n, k, i, j;
double alpha, beta;
printf ("\n This example computes real matrix C=alpha*A*B+beta*C using \n"
" Intel(R) MKL function dgemm, where A, B, and C are matrices and \n"
" alpha and beta are double precision scalars\n\n");
m = 2000, k = 200, n = 1000;
printf (" Initializing data for matrix multiplication C=A*B for matrix \n"
" A(%ix%i) and matrix B(%ix%i)\n\n", m, k, k, n);
alpha = 1.0; beta = 0.0;
printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n"
" performance \n\n");
A = (double *)malloc( m*k*sizeof( double ));
B = (double *)malloc( k*n*sizeof( double ));
C = (double *)malloc( m*n*sizeof( double ));
printf (" Intializing matrix data \n\n");
for (i = 0; i < (m*k); i++) {
A[i] = (double)(i+1);
}
for (i = 0; i < (k*n); i++) {
B[i] = (double)(-i-1);
}
for (i = 0; i < (m*n); i++) {
C[i] = 0.0;
}
printf (" Computing matrix product using Intel(R) MKL dgemm function via CBLAS interface \n\n");
// cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
// m, n, k, alpha, A, k, B, n, beta, C, n);
// eigenlib_dgemm(const char *TransA, const char *TransB, const double alpha,
// const gsl_matrix *A, const gsl_matrix *B, const double beta,
// gsl_matrix *C) {
gsl_matrix *AM = gsl_matrix_safe_alloc(m,k); // rows x cols
gsl_matrix *BM = gsl_matrix_safe_alloc(k,n);
gsl_matrix *CM = gsl_matrix_calloc(m,n);
fast_copy(AM,A);
fast_copy(BM,B);
fast_copy(CM,C);
fast_dgemm("N","N",alpha,AM,BM,beta,CM);
printf ("\n Computations completed.\n\n");
A = AM->data;
B = BM->data;
C = CM->data;
REQUIRE(trunc(C[0]) == -2666620100.0 );
REQUIRE(trunc(C[1]) == -2666640200.0 );
REQUIRE(trunc(C[2003]) == -10627000400.0 );
}
// The following code is normally disabled as it stops the unit tests!
#ifdef TEST_SIGFPE
#include <stdio.h>
#include <sys/types.h>
#include <unistd.h>
#include <signal.h>
void sighandler(int signum)
{
printf("Process %d got signal %d\n", getpid(), signum);
signal(signum, SIG_DFL);
kill(getpid(), signum);
}
TEST_CASE("NaN handlers", "[math]") {
signal(SIGFPE, sighandler);
feenableexcept(FE_INVALID |
FE_DIVBYZERO |
FE_OVERFLOW |
FE_UNDERFLOW);
double dirty = 0.0;
double nanval= 0.0/dirty;
printf("Succeeded! dirty=%lf, nanval=%lf\n",dirty,nanval);
}
#endif // TEST_SIGFPE
