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

Question

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.

Solution

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;
}

Comments

[Haatschii]

While the question is undoubtedly a bit vague, I too don't see why it should be off topic. The question for an example code snipped, performing a specific task (Matrix multiplication), using a specific tool (Blas/Lapack), is not likely "to attract opinionated answers and spam". I guess the downvotes are primary because one would expect the Internet is full of such simple examples and therefore on suspects the OP was to lazy to do some research before asking. However in case of Blas/Lapack examples, my experience is in fact that easy to understand examples are in fact hard to find.

[dreamer]

@Haatschii Thanks for the support. I really did try to research the topic before I asked but as you said, it's not as easy to find as it might seem. Thanks for your answer as well.

[dreamer]

Thank you! I will give it a try.

[dreamer]

Thanks for the extension of your answer!