C++ Integration

To get started calling any matrix-free function provided by primate, such sl_trace or lanczos, simply pass any type with a .shape() and .matvec() member functions defined like so:

class LinOp {
  int nrow, ncol;
  
  LinOp(int nr, int nc) : nrow(nr), ncol(nc) {}
  
  void matvec(const float* input, float* output) const {
    ... // implementation details 
  }

  void shape() const { return std::make_pair(nrow, ncol); }
}

It’s up to you to ensure shape() yields the correct size; primate will supply vectors to input of size .shape().second (number of columns) and guarantees the pointer to the output will be at least shape().first (number of rows), no more.

The LinearOperator concept

primate’s generic API is enabled through C++20 concepts. Thus, the more exact statement is that any type respecting the LinearOperator concept shown below can be passed:

using FP = std::floating_point; 
template < typename T, FP float_t = typename T::value_type >
concept LinearOperator = requires(T A, const float_t* v, float_t* o) {
  { A.matvec(v, o) }; // o = A v
  { A.shape() } -> std::convertible_to< std::pair< size_t, size_t > >;
};

An instance A of type T is said to support the LinearOperator concept if it has:

1. A method Av \mapsto o, with signature A.matvec(const float_t* v, float_t* o)
2. A method yielding (\mathrm{card}(o), \mathrm{card}(v)), with signatureA.shape() -> pair< ... >

shape() should yield a pair (n,m) representing the sizes of the output and input vectors, respectively. This corresponds to the number of rows and columns in the matrix setting.

Other Concepts

Depending on the problem at hand, the supplied operator may need to meet other constraints. Here’s a short list additional operator concepts:

Concept	Supports	Signature	Requires
LinearOperator	A v \mapsto o	A.matvec(v, o)	NA
AdjointOperator	A^T v \mapsto o	A.rmatvec(v, o)	LinearOperator
AdditiveOperator	o \gets o + \alpha Av	A.matvec_add(v, alpha, o)	LinearOperator
AffineOperator	Sets t s.t. A + tB	A.set_parameter(t)	LinearOperator

Respecting these constraints is opt-in: if your operator is symmetric and you only need access to the Lanczos method, then any satisfying the LinearOperator concept is sufficient.