Get Gentoo!
gentoo.org sites
gentoo.org Wiki Bugs Forums Packages
Planet Archives Sources
Infra Status
summaryrefslogtreecommitdiff
path: root/btl/accuracy
diff options
context:
space:
mode:
Diffstat (limited to 'btl/accuracy')
-rw-r--r--btl/accuracy/lapack/diff.hh47
-rw-r--r--btl/accuracy/lapack/lapack_LU.hh74
-rw-r--r--btl/accuracy/lapack/lapack_QR.hh74
-rw-r--r--btl/accuracy/lapack/lapack_STEV.hh52
-rw-r--r--btl/accuracy/lapack/lapack_SVD.hh57
-rw-r--r--btl/accuracy/lapack/lapack_SYEV.hh60
-rw-r--r--btl/accuracy/lapack/lapack_cholesky.hh54
-rw-r--r--btl/accuracy/lapack/main_lapack.cpp120
-rw-r--r--btl/accuracy/lapack/timer.hh67
-rw-r--r--btl/accuracy/main_blas.cpp253
-rw-r--r--btl/accuracy/sizes.hh59
11 files changed, 917 insertions, 0 deletions
diff --git a/btl/accuracy/lapack/diff.hh b/btl/accuracy/lapack/diff.hh
new file mode 100644
index 0000000..6679655
--- /dev/null
+++ b/btl/accuracy/lapack/diff.hh
@@ -0,0 +1,47 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef DIFF_HH
+#define DIFF_HH
+
+#include <vector>
+
+double diff(const int& N, const double* x, const double* test)
+{
+ std::vector<double> d(x, x+N);
+ const int iONE = 1;
+ const double alpha = -1.;
+ daxpy_(&N, &alpha, test, &iONE, &d[0], &iONE);
+ return dnrm2_(&N, &d[0], &iONE)/dnrm2_(&N, x, &iONE);
+}
+
+template<typename T>
+T diff(const std::vector<T> x, const std::vector<T> test)
+{
+ return diff(static_cast<const int&>(x.size()), &x[0], &test[0]);
+}
+
+//float diff(const int& N, const float* x, const float* test)
+//{
+// std::vector<float> d(x, x+N);
+// const int iONE = 1;
+// const float alpha = -1.;
+// saxpy_(&N, &alpha, test, &iONE, &d[0], &iONE);
+// return snrm2_(&N, &d[0], &iONE)/snrm2_(&N, x, &iONE);
+//}
+
+#endif /* DIFF_HH */
diff --git a/btl/accuracy/lapack/lapack_LU.hh b/btl/accuracy/lapack/lapack_LU.hh
new file mode 100644
index 0000000..686a97f
--- /dev/null
+++ b/btl/accuracy/lapack/lapack_LU.hh
@@ -0,0 +1,74 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef LAPACK_LU_HH
+#define LAPACK_LU_HH
+
+#include "LinearCongruential.hh"
+#include "diff.hh"
+
+vector<int> get_piv(const vector<int>& ipiv)
+{
+ const int size = ipiv.size();
+ vector<int> ret(size);
+ for (int i = 0; i < size; ++i)
+ ret[i] = i;
+
+ int ii, jj, tmp;
+ for (int i = 1; i <= size; ++i) {
+ ii = i-1;
+ jj = ipiv[ii]-1;
+ tmp = ret[ii];
+ ret[ii] = ret[jj];
+ ret[jj] = tmp;
+ }
+
+ return ret;
+}
+
+double test_LU(const int& N, const unsigned& seed = 0)
+{
+ LinearCongruential lc(seed);
+
+ vector<double> A(N*N), Acopy, P(N*N);
+
+ /* Fill A and Acopy */
+ for (int i = 0; i < N*N; ++i)
+ A[i] = lc.get_01();
+ Acopy = A;
+
+ /* Compute decomposition */
+ vector<int> ipiv(N);
+ int info;
+ dgetrf_(&N, &N, &A[0], &N, &ipiv[0], &info);
+ vector<int> piv = get_piv(ipiv);
+
+ /* Construct P */
+ for (int r = 0; r < N; ++r) {
+ int c = piv[r];
+ P[c+N*r] = 1;
+ }
+
+ /* Test */
+ const double alpha = 1.;
+ dtrmm_("R", "L", "N", "U", &N, &N, &alpha, &A[0], &N, &P[0], &N);
+ dtrmm_("R", "U", "N", "N", &N, &N, &alpha, &A[0], &N, &P[0], &N);
+
+ return diff(Acopy, P);
+}
+
+#endif /* LAPACK_LU_HH */
diff --git a/btl/accuracy/lapack/lapack_QR.hh b/btl/accuracy/lapack/lapack_QR.hh
new file mode 100644
index 0000000..8b9ba82
--- /dev/null
+++ b/btl/accuracy/lapack/lapack_QR.hh
@@ -0,0 +1,74 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef LAPACK_QR_HH
+#define LAPACK_QR_HH
+
+#include "LinearCongruential.hh"
+#include "diff.hh"
+
+double test_QR(const int& N, const unsigned& seed = 0)
+{
+ LinearCongruential lc(seed);
+ vector<double> A(N*N), Acopy, tau(N), work(1), Q(N*N), H(N*N), tmp;
+
+ /* Fill A and Acopy */
+ for (int i = 0; i < N*N; ++i)
+ A[i] = lc.get_01();
+ Acopy = A;
+
+ /* Retrieve lwork */
+ int lwork = -1, info;
+ dgeqrf_(&N, &N, &A[0], &N, &tau[0], &work[0], &lwork, &info);
+ lwork = work[0];
+ work.resize(lwork);
+
+ /* Compute decomposition */
+ dgeqrf_(&N, &N, &A[0], &N, &tau[0], &work[0], &lwork, &info);
+
+ /* Compute Q */
+ for (int i = 0; i < N; ++i)
+ Q[i+i*N] = 1.;
+
+ double alpha, t;
+ const double dZERO = 0., dONE = 1.;
+ const int iONE = 1, N2 = N*N;
+ work.resize(N2);
+ for (int i = 0; i < N; ++i) {
+ /* Generate H_i */
+ for (int r = 0; r < N; ++r)
+ for (int c = r; c < N; ++c)
+ H[r+c*N] = (r == c);
+ dcopy_(&N, &A[i*N], &iONE, &work[0], &iONE);
+ for (int j = 0; j < i; ++j)
+ work[j] = 0.;
+ work[i] = 1.;
+ t = -tau[i];
+ dsyr_("U", &N, &t, &work[0], &iONE, &H[0], &N);
+
+ /* Multiply Q = Q*H_i */
+ dsymm_("R", "U", &N, &N, &dONE, &H[0], &N, &Q[0], &N, &dZERO, &work[0], &N);
+ dcopy_(&N2, &work[0], &iONE, &Q[0], &iONE);
+ }
+
+ /* Multiply */
+ dtrmm_("R", "U", "N", "N", &N, &N, &dONE, &A[0], &N, &Q[0], &N);
+
+ return diff(Acopy, Q);
+}
+
+#endif /* LAPACK_QR_HH */
diff --git a/btl/accuracy/lapack/lapack_STEV.hh b/btl/accuracy/lapack/lapack_STEV.hh
new file mode 100644
index 0000000..3fc9147
--- /dev/null
+++ b/btl/accuracy/lapack/lapack_STEV.hh
@@ -0,0 +1,52 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef LAPACK_STEV_HH
+#define LAPACK_STEV_HH
+
+#include "LinearCongruential.hh"
+#include "diff.hh"
+
+double test_STEV(const int& N, const unsigned& seed = 0)
+{
+ LinearCongruential lc(seed);
+ vector<double> A(N*N), D(N), E(N-1), Z(N*N), DD(N*N), work(max(1, 2*N-2)), tmp(N*N);
+
+ /* Fill D, E and A */
+ for (int i = 0; i < N-1; ++i) {
+ D[i] = lc.get_01(); A[i+i*N] = D[i];
+ E[i] = lc.get_01(); A[(i+1)+i*N] = E[i]; A[i+(i+1)*N] = E[i];
+ }
+ D[N-1] = lc.get_01(); A[N*N-1] = D[N-1];
+
+ /* Perform computation */
+ int info;
+ dstev_("V", &N, &D[0], &E[0], &Z[0], &N, &work[0], &info);
+
+ /* Construct DD */
+ for (int i = 0; i < N; ++i)
+ DD[i+i*N] = D[i];
+
+ /* A = V*D*V' */
+ const double dONE = 1., dZERO = 0.;
+ dgemm_("N", "N", &N, &N, &N, &dONE, &Z[0], &N, &DD[0], &N, &dZERO, &tmp[0], &N);
+ dgemm_("N", "T", &N, &N, &N, &dONE, &tmp[0], &N, &Z[0], &N, &dZERO, &DD[0], &N);
+
+ return diff(A, DD);
+}
+
+#endif /* LAPACK_STEV_HH */
diff --git a/btl/accuracy/lapack/lapack_SVD.hh b/btl/accuracy/lapack/lapack_SVD.hh
new file mode 100644
index 0000000..7dc2ab1
--- /dev/null
+++ b/btl/accuracy/lapack/lapack_SVD.hh
@@ -0,0 +1,57 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef LAPACK_SVD_HH
+#define LAPACK_SVD_HH
+
+#include "LinearCongruential.hh"
+#include "diff.hh"
+
+double test_SVD(const int& N, const unsigned& seed = 0)
+{
+ LinearCongruential lc(seed);
+
+ vector<double> A(N*N), Acopy, U(N*N), VT(N*N), S_(N), S(N*N);
+ vector<double> work(1);
+
+ /* Fill A and Acopy */
+ for (int i = 0; i < N*N; ++i)
+ A[i] = lc.get_01();
+ Acopy = A;
+
+ /* Retrieve lwork */
+ int lwork = -1, info;
+ dgesvd_("A", "A", &N, &N, &A[0], &N, &S_[0], &U[0], &N, &VT[0], &N, &work[0], &lwork, &info);
+ lwork = work[0];
+ work.resize(lwork);
+
+ /* Compute decomposition */
+ dgesvd_("A", "A", &N, &N, &A[0], &N, &S_[0], &U[0], &N, &VT[0], &N, &work[0], &lwork, &info);
+
+ /* Construct S */
+ for (int r = 0; r < N; ++r)
+ S[r+r*N] = S_[r];
+
+ /* Test */
+ const double dONE = 1., dZERO = 0.;
+ dgemm_("N", "N", &N, &N, &N, &dONE, &U[0], &N, &S[0], &N, &dZERO, &A[0], &N);
+ dgemm_("N", "N", &N, &N, &N, &dONE, &A[0], &N, &VT[0], &N, &dZERO, &S[0], &N);
+
+ return diff(Acopy, S);
+}
+
+#endif /* LAPACK_SVD_HH */
diff --git a/btl/accuracy/lapack/lapack_SYEV.hh b/btl/accuracy/lapack/lapack_SYEV.hh
new file mode 100644
index 0000000..dc6b733
--- /dev/null
+++ b/btl/accuracy/lapack/lapack_SYEV.hh
@@ -0,0 +1,60 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef LAPACK_SYEV_HH
+#define LAPACK_SYEV_HH
+
+#include "LinearCongruential.hh"
+#include "diff.hh"
+
+double test_SYEV(const int& N, const unsigned& seed = 0)
+{
+ LinearCongruential lc(seed);
+ vector<double> A(N*N), Acopy, W(N), D(N*N), work(1), tmp(N*N);
+
+ /* Fill A (symmetric) and Acopy */
+ for (int r = 0; r < N; ++r) {
+ A[r+r*N] = lc.get_01();
+ for (int c = r+1; c < N; ++c) {
+ A[r+c*N] = lc.get_01();
+ A[c+r*N] = A[r+c*N];
+ }
+ }
+ Acopy = A;
+
+ /* Retrieve lwork */
+ int lwork = -1, info;
+ dsyev_("V", "U", &N, &A[0], &N, &W[0], &work[0], &lwork, &info);
+ lwork = work[0];
+ work.resize(lwork);
+
+ /* Perform computation */
+ dsyev_("V", "U", &N, &A[0], &N, &W[0], &work[0], &lwork, &info);
+
+ /* Construct D */
+ for (int i = 0; i < N; ++i)
+ D[i+i*N] = W[i];
+
+ /* A = V*D*V' */
+ const double dONE = 1., dZERO = 0.;
+ dgemm_("N", "N", &N, &N, &N, &dONE, &A[0], &N, &D[0], &N, &dZERO, &tmp[0], &N);
+ dgemm_("N", "T", &N, &N, &N, &dONE, &tmp[0], &N, &A[0], &N, &dZERO, &D[0], &N);
+
+ return diff(Acopy, D);
+}
+
+#endif /* LAPACK_SYEV_HH */
diff --git a/btl/accuracy/lapack/lapack_cholesky.hh b/btl/accuracy/lapack/lapack_cholesky.hh
new file mode 100644
index 0000000..c398c0a
--- /dev/null
+++ b/btl/accuracy/lapack/lapack_cholesky.hh
@@ -0,0 +1,54 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef LAPACK_CHOLESKY_HH
+#define LAPACK_CHOLESKY_HH
+
+#include "LinearCongruential.hh"
+#include "diff.hh"
+
+double test_cholesky(const int& N, const unsigned& seed = 0)
+{
+ LinearCongruential lc(seed);
+
+ vector<double> A(N*N), Acopy, test(N*N);
+
+ /* Fill A (SPD), Acopy and test (identity) */
+ for (int r = 0; r < N; ++r) {
+ A[r+r*N] = N;
+ test[r+r*N] = 1.;
+ for (int c = r; c < N; ++c) {
+ A[r+c*N] += lc.get_01();
+ A[c+r*N] = A[r+c*N];
+ }
+ }
+ Acopy = A;
+
+ /* Compute decomposition */
+ int info;
+ dpotrf_("L", &N, &A[0], &N, &info);
+
+ /* Test */
+ const double alpha = 1.;
+ dtrmm_("L", "L", "N", "N", &N, &N, &alpha, &A[0], &N, &test[0], &N);
+ dtrmm_("R", "L", "T", "N", &N, &N, &alpha, &A[0], &N, &test[0], &N);
+
+ return diff(Acopy, test);
+
+}
+
+#endif /* LAPACK_CHOLESKY_HH */
diff --git a/btl/accuracy/lapack/main_lapack.cpp b/btl/accuracy/lapack/main_lapack.cpp
new file mode 100644
index 0000000..bd5db21
--- /dev/null
+++ b/btl/accuracy/lapack/main_lapack.cpp
@@ -0,0 +1,120 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#include "lapack.hh"
+#include "LinearCongruential.hh"
+#include "timer.hh"
+#include "sizes.hh"
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <fstream>
+#include <sstream>
+#include <string>
+
+using namespace std;
+
+extern "C" {
+ void daxpy_(const int*, const double*, const double*, const int*, double*, const int*);
+ void dcopy_(const int*, const double*, const int*, double*, const int*);
+ double dnrm2_(const int*, const double*, const int*);
+ void dsyr_(const char*, const int*, const double*, const double*, const int*, double*, const int*);
+ void dger_(const int*, const int*, const double*, const double*, 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*);
+ void dsymm_(const char*, const char*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*);
+ void dtrmm_(const char*, const char*, const char*, const char*, const int*, const int*, const double*, const double*, const int*, double*, const int*);
+}
+#include "lapack_LU.hh"
+#include "lapack_cholesky.hh"
+#include "lapack_SVD.hh"
+#include "lapack_QR.hh"
+#include "lapack_SYEV.hh"
+#include "lapack_STEV.hh"
+
+template<typename exec_t>
+void test(exec_t exec, const std::string& testname, const int& max = 3000, const int& N = 100)
+{
+ vector<int> sizes = logsizes(1, max, N);
+ Timer timer;
+
+ ostringstream fname;
+ fname << "accuracy_" << testname << "_lapack.dat";
+ cout << "Testing " << testname << " --> " << fname.str() << endl;
+ ofstream fs(fname.str().c_str());
+
+ for (vector<int>::const_iterator i = sizes.begin(), end = sizes.end(); i != end; ++i) {
+ const int size = *i;
+ double error = 0;
+ timer.start();
+ int times = 0;
+ do
+ error += exec(size, times++);
+ while (timer.elapsed() < 1. || times < 16);
+ cout << " size: " << setw(4) << right << size;
+ cout << setw(15) << right << error/times << " ";
+ cout << "[" << setw(6) << times << " samples] ";
+ cout << "(" << setw(3) << right << (i-sizes.begin()+1) << "/" << N << ")" << endl;
+ fs << size << " " << error/times << "\n";
+ }
+
+ fs.close();
+}
+
+int main(int argc, char **argv)
+{
+ bool
+ lu_decomp=false, cholesky = false, svd_decomp=false, qr_decomp=false,
+ syev=false, stev=false
+ ;
+
+ for (int i = 1; i < argc; ++i) {
+ std::string arg = argv[i];
+ if (arg == "lu_decomp") lu_decomp = true;
+ else if (arg == "cholesky") cholesky = true;
+ else if (arg == "svd_decomp") svd_decomp = true;
+ else if (arg == "qr_decomp") qr_decomp = true;
+ else if (arg == "syev") syev = true;
+ else if (arg == "stev") stev = true;
+ }
+
+ if (lu_decomp) {
+ test(test_LU, "lu_decomp", 3000, 100);
+ }
+
+ if (cholesky) {
+ test(test_cholesky, "cholesky", 3000, 100);
+ }
+
+ if (svd_decomp) {
+ test(test_SVD, "svd_decomp", 1000, 90);
+ }
+
+ if (qr_decomp) {
+ test(test_QR, "qr_decomp", 500, 70);
+ }
+
+ if (syev) {
+ test(test_SYEV, "syev", 1000, 90);
+ }
+
+ if (stev) {
+ test(test_STEV, "stev", 1000, 90);
+ }
+
+ return 0;
+}
diff --git a/btl/accuracy/lapack/timer.hh b/btl/accuracy/lapack/timer.hh
new file mode 100644
index 0000000..82a0767
--- /dev/null
+++ b/btl/accuracy/lapack/timer.hh
@@ -0,0 +1,67 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef TIMER_HH
+#define TIMER_HH
+
+#include <sys/time.h>
+
+class Timer
+{
+public:
+ Timer() : ZERO(0.)
+ { }
+
+ void start()
+ {
+ t = walltime(&ZERO);
+ }
+
+ double elapsed()
+ {
+ return walltime(&t);
+ }
+
+private:
+ double t;
+
+ const double ZERO;
+
+ static double walltime(const double *t0)
+ {
+ double mic;
+ double time;
+ double mega = 0.000001;
+ struct timeval tp;
+ struct timezone tzp;
+ static long base_sec = 0;
+ static long base_usec = 0;
+
+ (void) gettimeofday(&tp, &tzp);
+ if (base_sec == 0) {
+ base_sec = tp.tv_sec;
+ base_usec = tp.tv_usec;
+ }
+
+ time = (double) (tp.tv_sec - base_sec);
+ mic = (double) (tp.tv_usec - base_usec);
+ time = (time + mic * mega) - *t0;
+ return (time);
+ }
+};
+
+#endif /* TIMER_HH */
diff --git a/btl/accuracy/main_blas.cpp b/btl/accuracy/main_blas.cpp
new file mode 100644
index 0000000..aceae58
--- /dev/null
+++ b/btl/accuracy/main_blas.cpp
@@ -0,0 +1,253 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#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;
+
+}
diff --git a/btl/accuracy/sizes.hh b/btl/accuracy/sizes.hh
new file mode 100644
index 0000000..be458bf
--- /dev/null
+++ b/btl/accuracy/sizes.hh
@@ -0,0 +1,59 @@
+//=====================================================
+// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+#ifndef SIZES_HH
+#define SIZES_HH
+
+#include <vector>
+#include <cmath>
+
+template<typename T>
+std::vector<T> logspace(const T& min, const T& max, const unsigned& N)
+{
+ std::vector<T> result;
+ result.reserve(N);
+
+ const double emin = std::log(static_cast<double>(min)),
+ emax = std::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>
+std::vector<T> logsizes(const T& min, const T& max, const unsigned& N)
+{
+ if (N <= 10)
+ return logspace(min, max, N);
+
+ std::vector<T> result;
+ result.reserve(N);
+
+ for (unsigned i = 0; i < 9; ++i)
+ result.push_back(i+1);
+ std::vector<T> lres = logspace(10, max, N-9);
+ for (unsigned i = 9; i < N; ++i)
+ result.push_back(lres[i-9]);
+
+ return result;
+}
+
+#endif /* SIZES_HH */