I’m physicist looking for a relatively ergonomic (and relatively efficient) way of computing the Pfaffian of a dense complex matrix in c++. For concreteness I would not be expecting to work with matrices larger than about 1000 by 1000 elements. I have attempted to use the package “pfafpack” however this is an extension of blas and the blas matrices are significantly less ergonomic to work with directly, for example I’d much prefer to write

C = A * B; 

than

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, k, B, n, beta, C, n);

Commonly used c++ libraries like boost and eigen do not appear to have this capability.