Matrix operations in C++ (using Blas/Lapack or some other alternative)

I am new to C++, and I am trying to learn how to do matrix operations in C++.

I have read that BLAS/LAPACK is the best way to do it (see http://cpplapack.sourceforge.net/). However, I find it difficult to get started with using it.

What is some example code on how I can do simple matrix operations, like matrix multiplication, inverses, etc., using BLAS/LAPACK in C++.

If it is easier using some other alternative, then I would also be curious to see some example code of that.

I assume that if you are new to C++, you are also new to C and Fortran. In that case I would definitely suggest to you, not to start with BLAS/LAPACK, at least not without a nice C++ wrapper.

My suggestion would be to have a look at Eigen which offers a much easier start to matrix operations using native C++ code. You can have a look at their tutorial to get started. The Eigen performance is said to be comparable to that of BLAS/LAPACK. See e.g. their benchmark. However I didn't test that myself.

If you really want to go low level and use BLAS/LAPACK, have a look at the available functions of cBlas (the C-Wrapper of BLAS) and LAPACK. Additionally, you can find some examples how to use Lapacke (The C-Wrapper of LAPACK) here. But don't expect things to be nice and well documented!

To finally give an answer to your question: Here is a code snipped I wrote some time ago for benchmarking. The code creates two random matrices A and B and multiplies them into the matrix C.

#include <random>
#include <cblas.h>

int main ( int argc, char* argv[] ) {

    // Random numbers
    std::mt19937_64 rnd;
    std::uniform_real_distribution<double> doubleDist(0, 1);

    // Create arrays that represent the matrices A,B,C
    const int n = 20;
    double*  A = new double[n*n];
    double*  B = new double[n*n];
    double*  C = new double[n*n];

    // Fill A and B with random numbers
    for(uint i =0; i <n; i++){
        for(uint j=0; j<n; j++){
            A[i*n+j] = doubleDist(rnd);
            B[i*n+j] = doubleDist(rnd);
        }
    }

    // Calculate A*B=C
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0, A, n, B, n, 0.0, C, n);

    // Clean up
    delete[] A;
    delete[] B;
    delete[] C;

    return 0;
}