Update / improvements

.gitignore

@@ -0,0 +1,6 @@
+*.swp
+*.pyc
+build/
+.lock-wafbuild
+.unittest-gtest/
+.waf-*/

README.md

@@ -12,8 +12,8 @@ If you want to install it from source code,
 gfortran, BLAS (ATLAS is better as its implementation) and LAPACK are
 required.
 
-arpaca_performance_test also depends on
-[pficommon](http://github.com/pfi/pficommon "pficommon").
+`arpaca_performance_test` and `arpaca_performance_plot` require a compiler with
+C++11 support.
 
 Installation
 ------------
@@ -22,6 +22,20 @@ Just copy arpaca.hpp where you like.
 Since the interface may be changed in the future,
 it is recommended to copy it to the local directory of your own project.
 
+To build the example programs and test suite:
+
+```bash
+$ ./waf configure
+$ ./waf
+# build and run tests
+$ ./waf --check
+# optionally install arpaca.hpp to /usr/local/include
+$ ./waf install
+```
+
+You will find the programs in the directory `build`. If `./waf configure` does
+not succeed try changing the include and lib paths in `wscript`.
+
 Typical Usage
 -------------
 

arpaca.hpp

@@ -4,7 +4,7 @@
 #include <algorithm>
 #include <stdexcept>
 #include <vector>
-#include <eigen3/Eigen/Core>
+#include <Eigen/Core>
 
 // ARPACK interface
 extern "C" {

arpaca_test.cpp

@@ -1,11 +1,9 @@
-#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
-
 #include <cstdlib>
 #include <algorithm>
 #include <iostream>
 #include <vector>
-#include <eigen3/Eigen/Dense>
-#include <eigen3/Eigen/Sparse>
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
 #include <gtest/gtest.h>
 #include "arpaca.hpp"
 

performance_main.cpp

@@ -1,124 +1,117 @@
-#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+/*
+ * Copyright (c) 2011 Seiya Tokui <[email protected]>
+ * Copyright (c) 2014 Burkhard Ritter <[email protected]>
+ * This code is distributed under the MIT license.
+ *
+ * Performance test for Arpaca's eigensolver.
+ *
+ * To compile this program: 
+ * g++ \
+ *    -std=c++11 \
+ *    -I [/path/to/eigen] \
+ *    -O3 \
+ *    -DNDEBUG \
+ *    performance_main.cpp \
+ *    -L [/path/to/libarpack.a] \
+ *    -larpack \
+ *    -o performance_main
+ *
+ * This assumes that arpaca.hpp is in the same directory.
+ */
 
-#include <cstdlib>
-#include <algorithm>
 #include <iostream>
-#include <stdexcept>
-#include <utility>
 #include <vector>
-#include <tr1/cstdint>
-#include <tr1/unordered_map>
-#include <eigen3/Eigen/Core>
-#include <eigen3/Eigen/Sparse>
-#include <pficommon/math/random.h>
-#include <pficommon/system/time_util.h>
+#include <random>
+#include <chrono>
+#include <Eigen/Core>
+#include <Eigen/Sparse>
 #include "arpaca.hpp"
 
-struct EqualFirst {
-  template<typename Pair>
-  bool operator()(const Pair& l, const Pair& r) const
-  {
-    return l.first == r.first;
-  }
-};
+namespace chrono = std::chrono;
+
+typedef Eigen::Triplet<double> T;
+typedef Eigen::SparseMatrix<double, Eigen::RowMajor> SMatrix;
 
 template<typename RandomNumberGenerator>
-Eigen::SparseMatrix<double, Eigen::RowMajor>
+SMatrix
 MakeSparseSymmetricRandomMatrix(int n, int k, RandomNumberGenerator& rnd)
 {
-  std::vector<std::pair<std::pair<int, int>, double> > elems;
-  elems.reserve(2 * n * k);
-
-  std::vector<int> seq;
-  for (int i = 0; i < n; ++i) {
-    pfi::math::random::sample_with_replacement(rnd, n, k, seq);
-    for (int j = 0; j < k; ++j)
-      elems.push_back(
-          std::make_pair(std::make_pair(i, seq[j]), rnd.next_gaussian()));
-  }
-
-  // make it symmetric
-  for (size_t i = 0, orig_size = elems.size(); i < orig_size; ++i) {
-    const std::pair<int, int>& pos = elems[i].first;
-    elems.push_back(
-        std::make_pair(std::make_pair(pos.second, pos.first), elems[i].second));
-  }
-
-  std::sort(elems.begin(), elems.end());
-  elems.erase(std::unique(elems.begin(), elems.end(), EqualFirst()),
-              elems.end());
-
-  Eigen::SparseMatrix<double, Eigen::RowMajor> mat(n, n);
-
-  int cur_row = 0;
-  for (size_t i = 0; i < elems.size(); ++i) {
-    const int row = elems[i].first.first;
-    const int col = elems[i].first.second;
-
-    if (cur_row < row) {
-      cur_row = row;
-      mat.startVec(row);
+    std::uniform_int_distribution<int> r_int(0,n-1);
+    std::normal_distribution<double> r_real;
+
+    std::vector<T> ts;
+    ts.reserve(n*k);
+    for (int l=0; l<n*k; l++)
+    {
+        // If we randomly create two or more triplets for the same matrix
+        // element, the values of all triplets will be summed...
+        int i = r_int(rnd);
+        int j = r_int(rnd);
+        double v = r_real(rnd);
+        ts.push_back(T(i,j,v));
     }
-    mat.insertBack(row, col) = elems[i].second;
-  }
 
-  mat.finalize();
-  return mat;
+    SMatrix mat(n, n);
+    mat.setFromTriplets(ts.begin(),ts.end());
+    return mat.selfadjointView<Eigen::Upper>();
 }
 
 arpaca::EigenvalueType
 GetEigenvalueType(const std::string& name)
 {
-  if (name == "LA")
-    return arpaca::ALGEBRAIC_LARGEST;
-  else if (name == "SA")
-    return arpaca::ALGEBRAIC_SMALLEST;
-  else if (name == "BE")
-    return arpaca::ALGEBRAIC_BOTH_END;
-  else if (name == "LM")
-    return arpaca::MAGNITUDE_LARGEST;
-  else if (name == "SM")
-    return arpaca::MAGNITUDE_SMALLEST;
-  throw std::invalid_argument("invalid eigenvalue type");
+    if (name == "LA")
+        return arpaca::ALGEBRAIC_LARGEST;
+    else if (name == "SA")
+        return arpaca::ALGEBRAIC_SMALLEST;
+    else if (name == "BE")
+        return arpaca::ALGEBRAIC_BOTH_END;
+    else if (name == "LM")
+        return arpaca::MAGNITUDE_LARGEST;
+    else if (name == "SM")
+        return arpaca::MAGNITUDE_SMALLEST;
+    throw std::invalid_argument("invalid eigenvalue type");
 }
 
 int main(int argc, char** argv)
 {
-  if (argc != 5) {
-    std::cerr << "usage: " << argv[0]
-              << " <dimension>"
-              << " <# of non-zero values in each row>"
-              << " <# of eigenvectors>"
-              << " <type of eigenvalues>"
-              << std::endl;
-    std::cerr << "\ttype of eigenvalues: LA SA BE LM SM" << std::endl;
-    return 1;
-  }
+    if (argc != 5) {
+        std::cerr << "usage: " << argv[0]
+                  << " <dimension>"
+                  << " <# of non-zero values in each row>"
+                  << " <# of eigenvectors>"
+                  << " <type of eigenvalues>"
+                  << std::endl;
+        std::cerr << "\ttype of eigenvalues: LA SA BE LM SM" << std::endl;
+        return 1;
+    }
+
+    const int n = std::atoi(argv[1]),
+    k = std::atoi(argv[2]),
+    r = std::atoi(argv[3]);
+    const arpaca::EigenvalueType type = GetEigenvalueType(argv[4]);
 
-  const int n = std::atoi(argv[1]),
-      k = std::atoi(argv[2]),
-      r = std::atoi(argv[3]);
-  const arpaca::EigenvalueType type = GetEigenvalueType(argv[4]);
+    std::cerr << "Making matrix" << std::endl;
+    std::mt19937 generator(42);
+    SMatrix X = MakeSparseSymmetricRandomMatrix(n, k, generator);
 
-  std::cerr << "Making matrix" << std::endl;
-  pfi::math::random::mtrand rnd;
-  Eigen::SparseMatrix<double, Eigen::RowMajor> X =
-      MakeSparseSymmetricRandomMatrix(n, k, rnd);
+    std::cerr << "Start performance test" << std::endl;
+    chrono::steady_clock::time_point begin = chrono::steady_clock::now();
 
-  std::cerr << "Start performance test" << std::endl;
-  pfi::system::time::clock_time begin = pfi::system::time::get_clock_time();
+    arpaca::SymmetricEigenSolver<double> solver = arpaca::Solve(X, r, type);
 
-  arpaca::SymmetricEigenSolver<double> solver = arpaca::Solve(X, r, type);
+    chrono::steady_clock::time_point end = chrono::steady_clock::now();
 
-  const double duration = pfi::system::time::get_clock_time() - begin;
+    chrono::duration<double> duration_ = 
+        chrono::duration_cast<chrono::duration<double>>(end - begin);
+    const double duration = duration_.count();
 
-  std::cout << "        DIMENSION: " << X.rows() << std::endl;
-  std::cout << "         NONZEROS: " << X.nonZeros() << std::endl;
-  std::cout << "         DURATION: "
-            << duration << " SEC." << std::endl;
-  std::cout << "             ITER: "
-            << solver.num_actual_iterations() << std::endl;
-  std::cout << "CONVERGED EIGVALS: "
-            << solver.num_converged_eigenvalues() << std::endl;
-  std::cout << "             INFO: " << solver.GetInfo() << std::endl;
+    std::cout << "        DIMENSION: " << X.rows() << std::endl;
+    std::cout << "         NONZEROS: " << X.nonZeros() << std::endl;
+    std::cout << "         DURATION: "
+              << duration << " SEC." << std::endl;
+    std::cout << "             ITER: "
+              << solver.num_actual_iterations() << std::endl;
+    std::cout << "CONVERGED EIGVALS: "
+              << solver.num_converged_eigenvalues() << std::endl;
+    std::cout << "             INFO: " << solver.GetInfo() << std::endl;
 }