diff options
author | 2012-10-21 15:52:07 +0200 | |
---|---|---|
committer | 2012-10-21 15:52:07 +0200 | |
commit | de6adbd75ee5ac479ae967df51fdf97991815042 (patch) | |
tree | 0a99c6b52a92db41d1716301d096b165bd16eb64 | |
parent | Added and tested module lapackAccuracy. (diff) | |
download | auto-numerical-bench-de6adbd75ee5ac479ae967df51fdf97991815042.tar.gz auto-numerical-bench-de6adbd75ee5ac479ae967df51fdf97991815042.tar.bz2 auto-numerical-bench-de6adbd75ee5ac479ae967df51fdf97991815042.zip |
Removed old accuracy framework.
-rw-r--r-- | btl/accuracy/lapack/diff.hh | 47 | ||||
-rw-r--r-- | btl/accuracy/lapack/lapack_LU.hh | 74 | ||||
-rw-r--r-- | btl/accuracy/lapack/lapack_QR.hh | 74 | ||||
-rw-r--r-- | btl/accuracy/lapack/lapack_STEV.hh | 52 | ||||
-rw-r--r-- | btl/accuracy/lapack/lapack_SVD.hh | 57 | ||||
-rw-r--r-- | btl/accuracy/lapack/lapack_SYEV.hh | 60 | ||||
-rw-r--r-- | btl/accuracy/lapack/lapack_cholesky.hh | 54 | ||||
-rw-r--r-- | btl/accuracy/lapack/main_lapack.cpp | 120 | ||||
-rw-r--r-- | btl/accuracy/lapack/timer.hh | 67 | ||||
-rw-r--r-- | btl/accuracy/main_blas.cpp | 253 | ||||
-rw-r--r-- | btl/accuracy/sizes.hh | 59 |
11 files changed, 0 insertions, 917 deletions
diff --git a/btl/accuracy/lapack/diff.hh b/btl/accuracy/lapack/diff.hh deleted file mode 100644 index 6679655..0000000 --- a/btl/accuracy/lapack/diff.hh +++ /dev/null @@ -1,47 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index 686a97f..0000000 --- a/btl/accuracy/lapack/lapack_LU.hh +++ /dev/null @@ -1,74 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index 8b9ba82..0000000 --- a/btl/accuracy/lapack/lapack_QR.hh +++ /dev/null @@ -1,74 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index 3fc9147..0000000 --- a/btl/accuracy/lapack/lapack_STEV.hh +++ /dev/null @@ -1,52 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index 7dc2ab1..0000000 --- a/btl/accuracy/lapack/lapack_SVD.hh +++ /dev/null @@ -1,57 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index dc6b733..0000000 --- a/btl/accuracy/lapack/lapack_SYEV.hh +++ /dev/null @@ -1,60 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index c398c0a..0000000 --- a/btl/accuracy/lapack/lapack_cholesky.hh +++ /dev/null @@ -1,54 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index bd5db21..0000000 --- a/btl/accuracy/lapack/main_lapack.cpp +++ /dev/null @@ -1,120 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index 82a0767..0000000 --- a/btl/accuracy/lapack/timer.hh +++ /dev/null @@ -1,67 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index aceae58..0000000 --- a/btl/accuracy/main_blas.cpp +++ /dev/null @@ -1,253 +0,0 @@ -//===================================================== -// 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 deleted file mode 100644 index be458bf..0000000 --- a/btl/accuracy/sizes.hh +++ /dev/null @@ -1,59 +0,0 @@ -//===================================================== -// 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 */ |