1 files changed, 236 insertions, 0 deletions
diff --git a/accuracy/main_blas.cpp b/accuracy/main_blas.cpp
new file mode 100644
index 0000000..13f8206
--- /dev/null
+++ b/accuracy/main_blas.cpp
@@ -0,0 +1,236 @@
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <cmath>
+#include <cstdlib>
+
+#define ADD_
+extern "C" {
+	void daxpy_(const int*, const double*, const double*, const int*, double*, const int*);
+	void dgemv_(const char*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*);
+	void dtrsv_(const char*, const char*, const char*, const int*, const double*, const int*, double*, const int*);
+	void dgemm_(const char*, const char*, const int*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*);
+	double dnrm2_(const int*, const double*, const int*);
+}
+
+using namespace std;
+
+template<typename T>
+void print_array(const int& N, const T* array) {
+	for (const T *p = array, *e = array+N; p != e; ++p)
+		cout << *p << " ";
+}
+
+template<typename T>
+void print_matrix(const int& rows, const int& cols, const T* matrix) {
+	for (int row = 0; row < rows; ++row) {
+		for (int col = 0; col < cols; ++col)
+			cout << *(matrix + rows*col + row) << " ";
+		cout << "\n";
+	}
+}
+
+template<typename T>
+vector<T> logspace(const T& min, const T& max, const unsigned& N)
+{
+	vector<T> result;
+	result.reserve(N);
+
+	const double emin = log(static_cast<double>(min)), emax = log(static_cast<double>(max));
+	double e, r = (emax-emin)/(N-1);
+	for (unsigned i = 0; i < N; ++i) {
+		e = emin + i*r;
+		result.push_back(static_cast<T>(exp(e)));
+	}
+
+	return result;
+}
+
+template<typename T>
+vector<T> logsizes(const T& min, const T& max, const unsigned& N)
+{
+	if (N <= 10)
+		return logspace(min, max, N);
+
+	vector<T> result;
+	result.reserve(N);
+
+	for (unsigned i = 0; i < 9; ++i)
+		result.push_back(i+1);
+	vector<T> lres = logspace(10, max, N-9);
+	for (unsigned i = 9; i < N; ++i)
+		result.push_back(lres[i-9]);
+
+	return result;
+}
+
+double axpy_test(const int& size)
+{
+	// Set up
+	const int ONE = 1;
+	const double alpha = 1./3.;
+	double *x = new double[size], *y = new double[size];
+	for (int i = 1; i <= size; ++i) {
+		x[i-1] = 10. / i;
+		y[i-1] = -10./(3. * i);
+	}
+
+	// Execute
+	daxpy_(&size, &alpha, x, &ONE, y, &ONE);
+
+	// Compute the error
+	double error = dnrm2_(&size, y, &ONE);
+
+	delete[] x; delete[] y;
+	return error;
+}
+
+double matrix_vector_test(const int& size)
+{
+	// Set up
+	const int ONE = 1;
+	char TRANS = 'N';
+	const double alpha = 1./size, beta = 1.;
+	double *A = new double[size*size], *x = new double[size], *y = new double[size];
+	for (int i = 1; i <= size; ++i) {
+		x[i-1] = i;
+		y[i-1] = -i;
+		for (int j = 1; j <= size; ++j) {
+			*(A+i-1+(j-1)*size) = static_cast<double>(i)/j;
+		}
+	}
+
+	// Execute
+	dgemv_(&TRANS, &size, &size, &alpha, A, &size, x, &ONE, &beta, y, &ONE);
+
+	// Compute the error
+	double error = dnrm2_(&size, y, &ONE);
+
+	delete[] A; delete[] x; delete[] y;
+	return error;
+}
+
+double trisolve_vector_test(const int& size)
+{
+	// Set up
+	const int ONE = 1;
+	char UPLO = 'U', TRANS = 'N', DIAG = 'U';
+	double *A = new double[size*size], *x = new double[size+1], *y = new double[size];
+	const double alpha = 1.;
+	x[size] = 0.;
+	for (int i = 1; i <= size; ++i) {
+		x[size-i] = x[size-i+1] + 1./i;
+		y[i-1] = -1.;
+		for (int j = i; j <= size; ++j) {
+			*(A+i-1+(j-1)*size) = 1. / (j-i+1);
+		}
+	}
+
+	// Execute
+	dtrsv_(&UPLO, &TRANS, &DIAG, &size, A, &size, x, &ONE);
+	daxpy_(&size, &alpha, x, &ONE, y, &ONE);
+	double error = dnrm2_(&size, y, &ONE);
+
+	delete[] A; delete[] x; delete[] y;
+	return error;
+}
+
+double matrix_matrix_test(const int& size)
+{
+	// rand48 initialization
+	srand48(5);
+
+	// sigma = SUM[k=1..size](k)
+	double sigma = 0;
+	for (int i = 1; i <= size; ++i)
+		sigma += i;
+
+	// Set up
+	const int ONE = 1;
+	char TRANS = 'N';
+	const double alpha = drand48(), beta = -2. * alpha * sigma;
+	const int size2 = size*size;
+	double *A = new double[size2], *B = new double[size2], *C = new double[size2];
+
+	for (int i = 1; i <= size; ++i)
+		for (int j = 1; j <= size; ++j) {
+			*(A+i-1+size*(j-1)) = static_cast<double>(j)/i;
+			*(B+i-1+size*(j-1)) = j;
+			*(C+i-1+size*(j-1)) = static_cast<double>(j)/(2.*i);
+		}
+
+	// Execute
+	dgemm_(&TRANS, &TRANS, &size, &size, &size, &alpha, A, &size, B, &size, &beta, C, &size);
+	double error = dnrm2_(&size2, C, &ONE);
+
+	delete[] A; delete[] B; delete[] C;
+	return error;
+}
+
+template<typename T>
+void test(T t, const int& min, const int& max, const unsigned& N, const string& name)
+{
+    ostringstream fname;
+	ofstream fout;
+	double result;
+	int N_;
+	
+	fname << "accuracy_" << name << "_blas.dat";
+	fout.open(fname.str().c_str());
+	cout << name << " test -- " << fname.str() << endl;
+	vector<int> axpy_sizes = logsizes(min, max, N);
+	N_ = 0;
+	for (vector<int>::const_reverse_iterator i = axpy_sizes.rbegin(), e = axpy_sizes.rend(); i!=e; ++i) {
+		result = t(*i);
+		fout << *i << " " << result << endl;
+		cout << " size: " << *i << "   " << result << "   (" << ++N_ << "/100)" << endl;
+	}
+	fout.close();
+	cout << "\n";
+}
+
+int main(int argv, char **argc)
+{
+    bool
+    axpy=false, axpby=false, rot=false,
+    matrix_vector=false, atv=false, symv=false, syr2=false, ger=false, trisolve_vector=false,
+    matrix_matrix=false, aat=false, trisolve_matrix=false, trmm=false  
+    ;
+  
+  
+    for (int i = 1; i < argv; ++i) {
+        std::string arg = argc[i];
+        if (arg == "axpy") axpy = true;
+        else if (arg == "axpby") axpby = true;
+        else if (arg == "rot") rot = true;
+        else if (arg == "matrix_vector") matrix_vector = true;
+        else if (arg == "atv") atv = true;
+        else if (arg == "symv") symv = true;
+        else if (arg == "syr2") syr2 = true;
+        else if (arg == "ger") ger = true;
+        else if (arg == "trisolve_vector") trisolve_vector = true;
+        else if (arg == "matrix_matrix") matrix_matrix = true;
+        else if (arg == "aat") aat = true;
+        else if (arg == "trisolve_matrix") trisolve_matrix = true;
+        else if (arg == "trmm") trmm = true;
+    }
+    
+    
+	/* AXPY test */
+	if (axpy)
+    test(axpy_test, 1, 1000000, 100, "axpy");
+    	
+    if (matrix_vector)
+    test(matrix_vector_test, 1, 3000, 100, "matrix_vector");
+	
+	if (trisolve_vector)
+    test(trisolve_vector_test, 1, 3000, 100, "trisolve_vector");
+    
+    if (matrix_matrix)
+    test(matrix_matrix_test, 1, 2000, 100, "matrix_matrix");
+
+    return 0;
+
+}