Usage examples: serial solvers

These examples show how to find a few smallest eigenvalues of a real symmetric matrix and their respective eigenvectors. The matrix is a banded one with the bandwidth equal to 5. Note that the matrix is never stored in memory and is defined only by the rule of how it acts on a vector.

Eigen 3

#include <cmath>
#include <iostream>

// This example shows how to use ezARPACK and the Eigen3 storage backend
// to partially diagonalize a large sparse symmetric matrix
// and find a number of its low-lying eigenvalues.

#include <ezarpack/arpack_solver.hpp>
#include <ezarpack/storages/eigen.hpp>
#include <ezarpack/version.hpp>

using namespace ezarpack;
using namespace Eigen;

// Size of the matrix
const int N = 10000;

// We are going to use a band matrix with this bandwidth
const int bandwidth = 5;

// The number of low-lying eigenvalues we want to compute
const int N_ev = 10;

int main() {

  // Print ezARPACK version
  std::cout << "Using ezARPACK version " << EZARPACK_VERSION << std::endl;

  // Construct a solver object for the symmetric case.
  // For the Eigen3 storage backend, other options would be
  // * `arpack_solver<ezarpack::Asymmetric, eigen_storage>' for general
  //   real matrices;
  // * `arpack_solver<ezarpack::Complex, eigen_storage>' for general
  //   complex matrices.
  using solver_t = arpack_solver<ezarpack::Symmetric, eigen_storage>;
  solver_t solver(N);

  // Specify parameters for the solver
  using params_t = solver_t::params_t;
  params_t params(N_ev,               // Number of low-lying eigenvalues
                  params_t::Smallest, // We want the smallest eigenvalues
                  true);              // Yes, we want the eigenvectors
                                      // (Ritz vectors) as well

  // Linear operator representing multiplication of a given vector by our matrix
  // The operator must act on the 'in' vector and store results in 'out'.
  auto matrix_op = [](solver_t::vector_const_view_t in,
                      solver_t::vector_view_t out) {
    out.fill(0); // Clear result

    // out_i = \sum_j A_{ij} in_j
    // A_{ij} = |i-j| / (1 + i + j), if |i-j| <= bandwidth, zero otherwise
    for(int i = 0; i < N; ++i) {
      int j_min = std::max(0, i - bandwidth);
      int j_max = std::min(N - 1, i + bandwidth);
      for(int j = j_min; j <= j_max; ++j) {
        out(i) += double(std::abs(i - j)) / (1 + i + j) * in(j);
      }
    }
  };

  // Run diagonalization!
  solver(matrix_op, params);

  // Number of converged eigenvalues
  std::cout << solver.nconv() << " out of " << params.n_eigenvalues
            << " eigenvalues have converged" << std::endl;

  // Print found eigenvalues
  std::cout << "Eigenvalues (Ritz values):" << std::endl;
  std::cout << solver.eigenvalues().transpose() << std::endl;

  // Check A*v = \lambda*v
  auto const& lambda = solver.eigenvalues();
  auto const& v = solver.eigenvectors();
  VectorXd lhs(N), rhs(N);

  for(int i = 0; i < N_ev; ++i) { // For each eigenpair ...
    const VectorXd eigenvec = v.col(i);
    matrix_op(eigenvec.head(N), lhs.head(N)); // calculate A*v
    rhs = lambda(i) * eigenvec;               // and \lambda*v

    std::cout << i << ": deviation = " << (rhs - lhs).norm() / N << std::endl;
  }

  // Print some computation statistics
  auto stats = solver.stats();

  std::cout << "Number of Lanczos update iterations: " << stats.n_iter
            << std::endl;
  std::cout << "Total number of OP*x operations: " << stats.n_op_x_operations
            << std::endl;
  std::cout << "Total number of steps of re-orthogonalization: "
            << stats.n_reorth_steps << std::endl;

  return 0;
}

Blaze

#include <cmath>
#include <iostream>

// This example shows how to use ezARPACK and the Blaze storage backend
// to partially diagonalize a large sparse symmetric matrix
// and find a number of its low-lying eigenvalues.

#include <ezarpack/arpack_solver.hpp>
#include <ezarpack/storages/blaze.hpp>
#include <ezarpack/version.hpp>

using namespace ezarpack;
using namespace blaze;

// Size of the matrix
const int N = 10000;

// We are going to use a band matrix with this bandwidth
const int bandwidth = 5;

// The number of low-lying eigenvalues we want to compute
const int N_ev = 10;

int main() {

  // Print ezARPACK version
  std::cout << "Using ezARPACK version " << EZARPACK_VERSION << std::endl;

  // Construct a solver object for the symmetric case.
  // For the Blaze storage backend, other options would be
  // * `arpack_solver<ezarpack::Asymmetric, blaze_storage>' for general
  //   real matrices;
  // * `arpack_solver<ezarpack::Complex, blaze_storage>' for general
  //   complex matrices.
  using solver_t = arpack_solver<ezarpack::Symmetric, blaze_storage>;
  solver_t solver(N);

  // Specify parameters for the solver
  using params_t = solver_t::params_t;
  params_t params(N_ev,               // Number of low-lying eigenvalues
                  params_t::Smallest, // We want the smallest eigenvalues
                  true);              // Yes, we want the eigenvectors
                                      // (Ritz vectors) as well

  // Linear operator representing multiplication of a given vector by our matrix
  // The operator must act on the 'in' vector and store results in 'out'.
  auto matrix_op = [](solver_t::vector_const_view_t in,
                      solver_t::vector_view_t out) {
    out.reset(); // Clear result

    // out_i = \sum_j A_{ij} in_j
    // A_{ij} = |i-j| / (1 + i + j), if |i-j| <= bandwidth, zero otherwise
    for(int i = 0; i < N; ++i) {
      int j_min = std::max(0, i - bandwidth);
      int j_max = std::min(N - 1, i + bandwidth);
      for(int j = j_min; j <= j_max; ++j) {
        out[i] += double(std::abs(i - j)) / (1 + i + j) * in[j];
      }
    }
  };

  // Run diagonalization!
  solver(matrix_op, params);

  // Number of converged eigenvalues
  std::cout << solver.nconv() << " out of " << params.n_eigenvalues
            << " eigenvalues have converged" << std::endl;

  // Print found eigenvalues
  std::cout << "Eigenvalues (Ritz values):" << std::endl;
  std::cout << trans(solver.eigenvalues()) << std::endl;

  // Check A*v = \lambda*v
  auto const& lambda = solver.eigenvalues();
  auto const& v = solver.eigenvectors();
  DynamicVector<double> lhs(N), rhs(N);

  for(int i = 0; i < N_ev; ++i) { // For each eigenpair ...
    const DynamicVector<double> eigenvec = column(v, i);
    matrix_op(subvector(eigenvec, 0, N), subvector(lhs, 0, N)); // calculate A*v
    rhs = lambda[i] * eigenvec;                                 // and \lambda*v

    std::cout << i << ": deviation = " << norm(rhs - lhs) / N << std::endl;
  }

  // Print some computation statistics
  auto stats = solver.stats();

  std::cout << "Number of Lanczos update iterations: " << stats.n_iter
            << std::endl;
  std::cout << "Total number of OP*x operations: " << stats.n_op_x_operations
            << std::endl;
  std::cout << "Total number of steps of re-orthogonalization: "
            << stats.n_reorth_steps << std::endl;

  return 0;
}

Armadillo

#include <cmath>
#include <iostream>

// This example shows how to use ezARPACK and the Armadillo storage backend
// to partially diagonalize a large sparse symmetric matrix
// and find a number of its low-lying eigenvalues.

#include <ezarpack/arpack_solver.hpp>
#include <ezarpack/storages/armadillo.hpp>
#include <ezarpack/version.hpp>

using namespace ezarpack;
using namespace arma;

// Size of the matrix
const int N = 10000;

// We are going to use a band matrix with this bandwidth
const int bandwidth = 5;

// The number of low-lying eigenvalues we want to compute
const int N_ev = 10;

int main() {

  // Print ezARPACK version
  std::cout << "Using ezARPACK version " << EZARPACK_VERSION << std::endl;

  // Construct a solver object for the symmetric case.
  // For the Armadillo storage backend, other options would be
  // * `arpack_solver<Asymmetric, armadillo_storage>' for general real matrices;
  // * `arpack_solver<Complex, armadillo_storage>' for general complex matrices.
  using solver_t = arpack_solver<Symmetric, armadillo_storage>;
  solver_t solver(N);

  // Specify parameters for the solver
  using params_t = solver_t::params_t;
  params_t params(N_ev,               // Number of low-lying eigenvalues
                  params_t::Smallest, // We want the smallest eigenvalues
                  true);              // Yes, we want the eigenvectors
                                      // (Ritz vectors) as well

  // Linear operator representing multiplication of a given vector by our matrix
  // The operator must act on the 'in' vector and store results in 'out'.
  auto matrix_op = [](solver_t::vector_const_view_t in,
                      solver_t::vector_view_t out) {
    out.zeros(); // Clear result

    // out_i = \sum_j A_{ij} in_j
    // A_{ij} = |i-j| / (1 + i + j), if |i-j| <= bandwidth, zero otherwise
    for(int i = 0; i < N; ++i) {
      int j_min = std::max(0, i - bandwidth);
      int j_max = std::min(N - 1, i + bandwidth);
      for(int j = j_min; j <= j_max; ++j) {
        out[i] += double(std::abs(i - j)) / (1 + i + j) * in[j];
      }
    }
  };

  // Run diagonalization!
  solver(matrix_op, params);

  // Number of converged eigenvalues
  std::cout << solver.nconv() << " out of " << params.n_eigenvalues
            << " eigenvalues have converged" << std::endl;

  // Print found eigenvalues
  std::cout << "Eigenvalues (Ritz values):" << std::endl;
  std::cout << solver.eigenvalues().t() << std::endl;

  // Check A*v = \lambda*v
  auto const& lambda = solver.eigenvalues();
  auto const& v = solver.eigenvectors();
  vec lhs(N), rhs(N);

  for(int i = 0; i < N_ev; ++i) { // For each eigenpair ...
    auto const eigenvec = v.col(i);
    matrix_op(eigenvec, lhs(span::all)); // calculate A*v
    rhs = lambda[i] * eigenvec;          // and \lambda*v

    std::cout << i << ": deviation = " << norm(rhs - lhs) / N << std::endl;
  }

  // Print some computation statistics
  auto stats = solver.stats();

  std::cout << "Number of Lanczos update iterations: " << stats.n_iter
            << std::endl;
  std::cout << "Total number of OP*x operations: " << stats.n_op_x_operations
            << std::endl;
  std::cout << "Total number of steps of re-orthogonalization: "
            << stats.n_reorth_steps << std::endl;

  return 0;
}

Boost uBLAS

#include <algorithm>
#include <cmath>
#include <iostream>

// This example shows how to use ezARPACK and the uBLAS storage backend
// to partially diagonalize a large sparse symmetric matrix
// and find a number of its low-lying eigenvalues.

#include <ezarpack/arpack_solver.hpp>
#include <ezarpack/storages/ublas.hpp>
#include <ezarpack/version.hpp>

#include <boost/numeric/ublas/io.hpp>

using namespace ezarpack;
using namespace boost::numeric::ublas;

// Size of the matrix
const int N = 10000;

// We are going to use a band matrix with this bandwidth
const int bandwidth = 5;

// The number of low-lying eigenvalues we want to compute
const int N_ev = 10;

int main() {

  // Print ezARPACK version
  std::cout << "Using ezARPACK version " << EZARPACK_VERSION << std::endl;

  // Construct a solver object for the symmetric case.
  // For the uBLAS storage backend, other options would be
  // * `arpack_solver<Asymmetric, ublas_storage>' for general real matrices;
  // * `arpack_solver<Complex, ublas_storage>' for general complex matrices.
  using solver_t = arpack_solver<Symmetric, ublas_storage>;
  solver_t solver(N);

  // Specify parameters for the solver
  using params_t = solver_t::params_t;
  params_t params(N_ev,               // Number of low-lying eigenvalues
                  params_t::Smallest, // We want the smallest eigenvalues
                  true);              // Yes, we want the eigenvectors
                                      // (Ritz vectors) as well

  // Linear operator representing multiplication of a given vector by our matrix
  // The operator must act on the 'in' vector and store results in 'out'.
  auto matrix_op = [](solver_t::vector_const_view_t in,
                      solver_t::vector_view_t out) {
    std::fill(out.begin(), out.end(), .0); // Clear result

    // out_i = \sum_j A_{ij} in_j
    // A_{ij} = |i-j| / (1 + i + j), if |i-j| <= bandwidth, zero otherwise
    for(int i = 0; i < N; ++i) {
      int j_min = std::max(0, i - bandwidth);
      int j_max = std::min(N - 1, i + bandwidth);
      for(int j = j_min; j <= j_max; ++j) {
        out(i) += double(std::abs(i - j)) / (1 + i + j) * in(j);
      }
    }
  };

  // Run diagonalization!
  solver(matrix_op, params);

  // Number of converged eigenvalues
  std::cout << solver.nconv() << " out of " << params.n_eigenvalues
            << " eigenvalues have converged" << std::endl;

  // Print found eigenvalues
  std::cout << "Eigenvalues (Ritz values):" << std::endl;
  std::cout << solver.eigenvalues() << std::endl;

  // Check A*v = \lambda*v
  auto const& lambda = solver.eigenvalues();
  auto const& v = solver.eigenvectors();
  vector<double> lhs(N), rhs(N);

  for(int i = 0; i < N_ev; ++i) { // For each eigenpair ...
    const vector<double> eigenvec = column(v, i);
    matrix_op(subrange(eigenvec, 0, N), subrange(lhs, 0, N)); // calculate A*v
    rhs = lambda(i) * eigenvec;                               // and \lambda*v

    std::cout << i << ": deviation = " << norm_2(rhs - lhs) / N << std::endl;
  }

  // Print some computation statistics
  auto stats = solver.stats();

  std::cout << "Number of Lanczos update iterations: " << stats.n_iter
            << std::endl;
  std::cout << "Total number of OP*x operations: " << stats.n_op_x_operations
            << std::endl;
  std::cout << "Total number of steps of re-orthogonalization: "
            << stats.n_reorth_steps << std::endl;

  return 0;
}

TRIQS arrays

#include <cmath>
#include <iostream>

// This example shows how to use ezARPACK and the TRIQS storage backend
// to partially diagonalize a large sparse symmetric matrix
// and find a number of its low-lying eigenvalues.

#include <ezarpack/arpack_solver.hpp>
#include <ezarpack/storages/triqs.hpp>
#include <ezarpack/version.hpp>

using namespace ezarpack;
using namespace triqs::arrays;

// Size of the matrix
const int N = 10000;

// We are going to use a band matrix with this bandwidth
const int bandwidth = 5;

// The number of low-lying eigenvalues we want to compute
const int N_ev = 10;

int main() {

  // Print ezARPACK version
  std::cout << "Using ezARPACK version " << EZARPACK_VERSION << std::endl;

  // Construct a solver object for the symmetric case.
  // For the TRIQS storage backend, other options would be
  // * `arpack_solver<Asymmetric, triqs_storage>' for general real matrices;
  // * `arpack_solver<Complex, triqs_storage>' for general complex matrices.
  using solver_t = arpack_solver<Symmetric, triqs_storage>;
  solver_t solver(N);

  // Specify parameters for the solver
  using params_t = solver_t::params_t;
  params_t params(N_ev,               // Number of low-lying eigenvalues
                  params_t::Smallest, // We want the smallest eigenvalues
                  true);              // Yes, we want the eigenvectors
                                      // (Ritz vectors) as well

  // Linear operator representing multiplication of a given vector by our matrix
  // The operator must act on the 'in' vector and store results in 'out'.
  auto matrix_op = [](auto in, auto out) {
    out() = 0; // Clear result

    // out_i = \sum_j A_{ij} in_j
    // A_{ij} = |i-j| / (1 + i + j), if |i-j| <= bandwidth, zero otherwise
    for(int i = 0; i < N; ++i) {
      int j_min = std::max(0, i - bandwidth);
      int j_max = std::min(N - 1, i + bandwidth);
      for(int j = j_min; j <= j_max; ++j) {
        out(i) += double(std::abs(i - j)) / (1 + i + j) * in(j);
      }
    }
  };

  // Run diagonalization!
  solver(matrix_op, params);

  // Number of converged eigenvalues
  std::cout << solver.nconv() << " out of " << params.n_eigenvalues
            << " eigenvalues have converged" << std::endl;

  // Print found eigenvalues
  std::cout << "Eigenvalues (Ritz values):" << std::endl;
  std::cout << solver.eigenvalues() << std::endl;

  // Check A*v = \lambda*v
  auto const& lambda = solver.eigenvalues();
  auto const& v = solver.eigenvectors();
  vector<double> lhs(N), rhs(N);

  for(int i = 0; i < N_ev; ++i) { // For each eigenpair ...
    auto const eigenvec = v(range(), i);
    matrix_op(eigenvec, lhs()); // calculate A*v
    rhs = lambda(i) * eigenvec; // and \lambda*v

    std::cout << i << ": deviation = " << norm2(rhs - lhs) / N << std::endl;
  }

  // Print some computation statistics
  auto stats = solver.stats();

  std::cout << "Number of Lanczos update iterations: " << stats.n_iter
            << std::endl;
  std::cout << "Total number of OP*x operations: " << stats.n_op_x_operations
            << std::endl;
  std::cout << "Total number of steps of re-orthogonalization: "
            << stats.n_reorth_steps << std::endl;

  return 0;
}

TRIQS/nda

#include <cmath>
#include <iostream>

// This example shows how to use ezARPACK and the TRIQS/nda storage backend
// to partially diagonalize a large sparse symmetric matrix
// and find a number of its low-lying eigenvalues.

#include <ezarpack/arpack_solver.hpp>
#include <ezarpack/storages/nda.hpp>
#include <ezarpack/version.hpp>

using namespace ezarpack;
using namespace nda;

// Size of the matrix
const int N = 10000;

// We are going to use a band matrix with this bandwidth
const int bandwidth = 5;

// The number of low-lying eigenvalues we want to compute
const int N_ev = 10;

int main() {

  // Print ezARPACK version
  std::cout << "Using ezARPACK version " << EZARPACK_VERSION << std::endl;

  // Construct a solver object for the symmetric case.
  // For the nda storage backend, other options would be
  // * `arpack_solver<Asymmetric, nda_storage>' for general real matrices;
  // * `arpack_solver<Complex, nda_storage>' for general complex matrices.
  using solver_t = arpack_solver<Symmetric, nda_storage>;
  solver_t solver(N);

  // Specify parameters for the solver
  using params_t = solver_t::params_t;
  params_t params(N_ev,               // Number of low-lying eigenvalues
                  params_t::Smallest, // We want the smallest eigenvalues
                  true);              // Yes, we want the eigenvectors
                                      // (Ritz vectors) as well

  // Linear operator representing multiplication of a given vector by our matrix
  // The operator must act on the 'in' vector and store results in 'out'.
  auto matrix_op = [](auto in, auto out) {
    out() = 0; // Clear result

    // out_i = \sum_j A_{ij} in_j
    // A_{ij} = |i-j| / (1 + i + j), if |i-j| <= bandwidth, zero otherwise
    for(int i = 0; i < N; ++i) {
      int j_min = std::max(0, i - bandwidth);
      int j_max = std::min(N - 1, i + bandwidth);
      for(int j = j_min; j <= j_max; ++j) {
        out(i) += double(std::abs(i - j)) / (1 + i + j) * in(j);
      }
    }
  };

  // Run diagonalization!
  solver(matrix_op, params);

  // Number of converged eigenvalues
  std::cout << solver.nconv() << " out of " << params.n_eigenvalues
            << " eigenvalues have converged" << std::endl;

  // Print found eigenvalues
  std::cout << "Eigenvalues (Ritz values):" << std::endl;
  std::cout << solver.eigenvalues() << std::endl;

  // Check A*v = \lambda*v
  auto const& lambda = solver.eigenvalues();
  auto const& v = solver.eigenvectors();
  vector<double> lhs(N), rhs(N);

  for(int i = 0; i < N_ev; ++i) { // For each eigenpair ...
    auto const eigenvec = v(range::all, i);
    matrix_op(eigenvec, lhs()); // calculate A*v
    rhs = lambda(i) * eigenvec; // and \lambda*v

    std::cout << i << ": deviation = " << std::sqrt(sum(abs2(rhs - lhs))) / N
              << std::endl;
  }

  // Print some computation statistics
  auto stats = solver.stats();

  std::cout << "Number of Lanczos update iterations: " << stats.n_iter
            << std::endl;
  std::cout << "Total number of OP*x operations: " << stats.n_op_x_operations
            << std::endl;
  std::cout << "Total number of steps of re-orthogonalization: "
            << stats.n_reorth_steps << std::endl;

  return 0;
}

xtensor

#include <cmath>
#include <iostream>

// This example shows how to use ezARPACK and the xtensor storage backend
// to partially diagonalize a large sparse symmetric matrix
// and find a number of its low-lying eigenvalues.

#include <ezarpack/arpack_solver.hpp>
#include <ezarpack/storages/xtensor.hpp>
#include <ezarpack/version.hpp>

#include <xtensor/io/xio.hpp>
#include <xtensor/reducers/xnorm.hpp>

using namespace ezarpack;
using namespace xt;

// Size of the matrix
const int N = 10000;

// We are going to use a band matrix with this bandwidth
const int bandwidth = 5;

// The number of low-lying eigenvalues we want to compute
const int N_ev = 10;

int main() {

  // Print ezARPACK version
  std::cout << "Using ezARPACK version " << EZARPACK_VERSION << std::endl;

  // Construct a solver object for the symmetric case.
  // For the xtensor storage backend, other options would be
  // * `arpack_solver<ezarpack::Asymmetric, xtensor_storage>' for general
  //   real matrices;
  // * `arpack_solver<ezarpack::Complex, xtensor_storage>' for general
  //   complex matrices.
  using solver_t = arpack_solver<ezarpack::Symmetric, xtensor_storage>;
  solver_t solver(N);

  // Specify parameters for the solver
  using params_t = solver_t::params_t;
  params_t params(N_ev,               // Number of low-lying eigenvalues
                  params_t::Smallest, // We want the smallest eigenvalues
                  true);              // Yes, we want the eigenvectors
                                      // (Ritz vectors) as well

  // Linear operator representing multiplication of a given vector by our matrix
  // The operator must act on the 'in' vector and store results in 'out'.
  auto matrix_op = [](auto in, auto out) {
    out.fill(.0); // Clear result

    // out_i = \sum_j A_{ij} in_j
    // A_{ij} = |i-j| / (1 + i + j), if |i-j| <= bandwidth, zero otherwise
    for(int i = 0; i < N; ++i) {
      int j_min = std::max(0, i - bandwidth);
      int j_max = std::min(N - 1, i + bandwidth);
      for(int j = j_min; j <= j_max; ++j) {
        out(i) += double(std::abs(i - j)) / (1 + i + j) * in(j);
      }
    }
  };

  // Run diagonalization!
  solver(matrix_op, params);

  // Number of converged eigenvalues
  std::cout << solver.nconv() << " out of " << params.n_eigenvalues
            << " eigenvalues have converged" << std::endl;

  // Print found eigenvalues
  std::cout << "Eigenvalues (Ritz values):" << std::endl;
  std::cout << solver.eigenvalues() << std::endl;

  // Check A*v = \lambda*v
  auto const& lambda = solver.eigenvalues();
  auto const& v = solver.eigenvectors();

  auto lhs = xt::xtensor<double, 1>::from_shape({N});
  auto rhs = xt::xtensor<double, 1>::from_shape({N});

  for(int i = 0; i < N_ev; ++i) { // For each eigenpair ...
    auto const eigenvec = xt::view(v, xt::all(), i);
    matrix_op(eigenvec, xt::view(lhs, xt::all())); // calculate A*v
    rhs = lambda(i) * eigenvec;                    // and \lambda*v

    std::cout << i << ": deviation = " << xt::norm_l2(rhs - lhs) / N
              << std::endl;
  }

  // Print some computation statistics
  auto stats = solver.stats();

  std::cout << "Number of Lanczos update iterations: " << stats.n_iter
            << std::endl;
  std::cout << "Total number of OP*x operations: " << stats.n_op_x_operations
            << std::endl;
  std::cout << "Total number of steps of re-orthogonalization: "
            << stats.n_reorth_steps << std::endl;

  return 0;
}