Generate a template class for a Matrix using vectors or change to Eigen library, C++

[baezacaljo]

Hello dear all i am new using this forum.

i have a question look, i have a project on which we are using several classes, therefore also several "Matrices", those matrices are composed using a vectors of vectors. Then, we perform several arithmetic operation like inverse multiplication and some additions. Hence, we decided to use Eigen library for such calculations, and in here is were the question appears (let tell you that by now everything works quite nice)

We are planning to keep doing some many more calculations because the project will grow, in other words we are going to generate some more classes that also require matrix operations.

Is the <vector> library the most suitable for such cases, for the generation of Matrices?. (BTW the size of the matrices are quite big some times goes from 10000 X 10000 up to 100000). if your answer is yes to keep using vectors, then my next question arises.

Do we have to generate a template class for the vector of vectors operations? or do we have to change to Eigen (sometimes eigen does not give zero values and gives small values close to zero). I've searched in the web and there are some guys that already tried to generate the matrix "solver" but they have huge problems with the efficiency and optimization. i know that we are far from perfect and in the it can be improved, but i am looking for your opinion experts.

However, if you say that we have to generate our own code can you address me through the right path, just give some ideas or things to read.

Thanks in advance.

[thmm]

I am not sure if you have to reinvent the wheel.
Have a look at the blaze library - it might have what you need.
https://bitbucket.org/blaze-lib/blaze/src/master/
https://www.youtube.com/watch?v=w-Y22KrMgFE

[lastchance]

(BTW the size of the matrices are quite big some times goes from 10000 X 10000 up to 100000)

You don't say what the applications of your matrices are, but there are very few where matrices of that size wouldn't be sparse (have only a limited number of non-zero elements). These tend to have specialist techniques for inversion.

You probably need access to sparse-matrix systems. Unless you have special types of matrix (e.g. tridiagonal or pentadiagonal for structured meshes) you would probably be better off sticking with Eigen.

[JLBorges]

> Is the <vector> library the most suitable for such cases, for the generation of Matrices?.

Continue using Eigen. Eigen's SparseMatrix if it is sparse: https://eigen.tuxfamily.org/dox/group__SparseQuickRefPage.html

[kigar64551]

Have you considered using MKL?
https://en.wikipedia.org/wiki/Math_Kernel_Library


See also:

https://www.intel.com/content/www/us/en/develop/documentation/mkl-tutorial-c/top/multiplying-matrices-using-dgemm.html

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/lapack-routines/lapack-linear-equation-routines/lapack-linear-equation-computational-routines/matrix-inversion-lapack-computational-routines.html

[jonnin]

eigen uses raw pointers, so if you need to do so you can do your own thing off to the side and use eigen when you need to as well, mixing and matching. There is probably a way to avoid copying the data when you do this, but I have not experimented.
do you actually NEED to be able to directly mess with the data outside of the objects eigen provides?

[baezacaljo]

i am just taking a look at blaze and i did not know anything about MKL.

About the sparse matrix, at this moment i have not used a matrix with most the elements containing zero's.

I am using eigen basically to perform the inverse and multiplication, that's the reason of my question i am looking for your opinion most of you are experts in C++, and know when something is already well optimized

[jonnin]

unless its special purpose any of these is probably fine.

most of these packages have some smarts and pick from several multiplies for the best one, saving you a lot of aggravation.

the thing about matrix math is that if you know one or two things about a matrix, you can do things faster that you can't do if you don't know them. An example is simply the eigenvalues... knowing those, you can avoid a lot of iteration for a number of algorithms that have to crunch hard without them. But getting them costs as much as what you were asking for. Catch 22... if you have them from earlier, great, but getting them saves nothing until you make use of them at least twice.

I don't know how these packages hand that sort of thing. Maybe they retain reusable shortcut info if you compute it and maybe they do not. Someone who uses this stuff more than I do may know.
But that is 'deep end of the pool' stuff too. If you are not deep in some project where every nanosecond counts, just let the library do what it does and be happy.

[kigar64551]

Whatever you do, a specialized library for matrix operations will probably be many times faster than a "naïve" implementation. MKL is well-known and widely-used, and it is heavily optimized for x86/x64 processors. But you can use whichever library gets the job done for you. Personally, I haven't used Eigen, only MKL.

The crux is: Most libraries that deal with matrices require those matrices to be stored in a specific way. The libraries often provide specialized data types and functions for that purpose. There can even be different data types for different "types" of matrices. For example, "sparse" matrices are very often stored in a CRS (Compressed Row Storage) or CSR (Compressed Sparse Row) format. Since it would be really bad if you had to convert, back and forth, between "your" matrix representation and the matrix representation used/required by the library all the time, I think you should model "your" matrix representation in such a way that it is "compatible" with the library you intend to use. Or just use the matrix data types from the library right away!

For example, the Eigen library apparently has a ready-to-use Matrix class. Why not use it then?
https://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html

those matrices are composed using a vectors of vectors

Now that would be a very inefficient way to store matrices, I think.

You would require one separate std::vector object for each row, and then another std::vector object where all those row-vectors are stored. Or you would require one separate std::vector object for each column, and then another std::vector object where all those column-vectors are stored. Either way, you get a whole lot of separate std::vector objects with a lot of separate "small" memory allocations and unnecessary overhead.

Furthermore: Almost certainly, you'd have to "convert" your matrix in order to pass it into the library.


It is probably better to use just a single std::vector object to store the entire matrix in one continuous chunk of memory. Individual elements can be addressed using either row-major indexing or column-major indexing:

https://en.wikipedia.org/wiki/Row-_and_column-major_order

But even that would only make sense for "dense" matrices, I suppose.

As mentioned before, for "sparse" matrices, special optimized formats, such as CRS/CSR, are usually used.

[baezacaljo]

So, with all that info, is much better to go ahead and do the magic with the library of my choice could be eigen MKL or even blaze. as far as i can see is better to use the library since the begging of the project right? because changing data from one vector of vectors to one Matrix of eigen library or MKL will take memory and of course time, it means copying and that is quite inefficient.

Also, as Jonnin mentioned, i am not a picky guy who is taking a look a millisecond difference, my project is not about that, even i can wait a bit longer for my result to come out of the solver. Also ill use big matrices and those libraries will provide the results

Finally i am using linux, and i am quite interested in MKL can you share some info? i just want to installed it and some tutorials.

Thanks a lot guys, now everything is clear.

[kigar64551]

as far as i can see is better to use the library since the begging of the project right? because changing data from one vector of vectors to one Matrix of eigen library or MKL will take memory and of course time, it means copying and that is quite inefficient.

That was my point, yes.

(plus: storing a matrix as a std::vector of many std::vector's is very inefficient, no matter what)


Best, decide for a library that you want to use, then use that libraries' data types (classes) all the way. Either that, or build a thin wrapper around the libraries' data types (classes). Avoid copying data back and forth.

i am quite interested in MKL can you share some info?

Intel has tons of documentation...

https://www.intel.com/content/www/us/en/develop/documentation/get-started-with-mkl-for-dpcpp/top.html

[jonnin]

when I used MKL it was just an optimized thin wrapper for the old stuff, Lapack and whatever... not sure what it is now. I didn't like those routines, you had to run like 5 pre-processing things to do a multiply.

[baezacaljo]

Tank you very much Kigar64551 and Jonnin.

I already started to do some modifications in the my code.

Now, i have a question, may be is quite silly and for beginners, but, what is a thin wrapper? and how to build it?.

may be that can lower my amount of work.

thanks for your answers and opinions.

[kigar64551]

Now, i have a question, may be is quite silly and for beginners, but, what is a thin wrapper? and how to build it?.

Suppose you decide to go with the Eigen library.

Now, you could just use Eigen's Matrix class everywhere in your code. But that would make your code base very dependent on Eigen library. Switching to another library later becomes hard.

Instead, you could create your own class MyMatrix. Internally "your" class would not re-invent the wheel but use an Eigen Matrix object to store everything; a "thin" layer of your own code around the "3rd-party" class.

The advantage is that your dependencies to the Eigen library are now "encapsulated" in your MyMatrix class. In case you ever decide you want to use a different library, changes only need to be done in one place.


It is even possible to create an interface, e.g. IMyMatrix, and then have several implementations (sub-classes) of that interface, such as MyEigenMatrix, MyMklMatrix and so on, for different libraries...

[jonnin]

what is a thin wrapper
a slight change to existing code. Usually you just change the interface, for example get rid of 3 or 4 parameters in a call to a function by creating a new function
void newfoo()
{
oldfoo(a,b,c); //where a,b,c are class variables in your wrapper now, presumably with appropriate values etc.
}

sometimes you add a tiny bit of functionality, but mostly it just reshaping stuff to make it easier to use/more consistent etc. MKL specifically rewrote or rebuilt a lot of it in assembly for intel chips, again, if I remember that part right, making it a faster version than the generic source code one.

[baezacaljo]

thanks kigar64551 and jonnin.

I got the whole point, basically is an small code that "links" the Eigen library and the classes of my project. As you said, if in the future i want to use another library only the "linker" "thin wrapper" is modified, which reduces the amount of work.

With that being said, i am fully interested in such thin wrapper class code. Therefore, i was searching on the web, and i can not get any nice and clear idea.

I came out with the idea of doing the thin wrapper using templates, on which i have to pass the the data

#pragma once
#include <bits/stdc++.h>
#include <Eigen/LU>
#include <string>

using namespace std;

template <typename T, const unsigned int _nRows, const unsigned int _nCols, const string _major>
class MATRIX
{
private:
    //The use of Eigen library, T is the typename it could be int or double or even float
    // _nRows is the number of rows in the matrix
    // _nCols is the number of columns in the matrix
    //_major is for the use of "row" or "column" major (how the matrix is place in the memory)
    Matrix<T, _nRows, _nCols, _major> MatrixWrapper;
public:
    // Constructor  
    MATRIX();
    // Copy constructor
    MATRIX(const MATRIX<T, _nRows, _nCols, _major>& other);
    // Destructor
    ~MATRIX();
    /*-------------- MATRIX/MATRIX Arithmetic --------------*/
    MATRIX<T, _nRows, _nCols, _major> operator*(const MATRIX<T, _nRows, _nCols, _major>& other);
    MATRIX<T, _nRows, _nCols, _major>& operator=(const MATRIX<T, _nRows, _nCols, _major>& other);

    /*-------------- Misc. MATRIX Operations --------------*/

    // Allows a matrix to be printed as output
    template<typename U, const unsigned int _nRowsPrint, const unsigned int _nColsPrint, const string _majorPrint>
    friend ostream& operator<<(ostream& print, const MATRIX<U, _nRowsPrint, _nColsPrint, _majorPrint>& MatrixWrapperPrint);

    /*-------------- Accessors --------------*/
    
    // Access to individual elements
    T& operator()(const unsigned int& nRow, const unsigned int& nCol);

    // Access to individual elements (constant)
    const T& operator()(const unsigned int& nRow, const unsigned int& nCol) const;

    // gets the number of rows
    unsigned int NumRows() const;

    // gets the number ofcolumns
    unsigned int NumCols() const;
};

Look, i know there is not a minimal reproducible example, but am i on track? or can you provide a minimal example?

Also, i am not use to work with using namespace std; and using namespace Eigen; just for a quick example is just wrote them.

Thank you very much guys you are helping me a lot

[jonnin]

I was talking about the MKL not eigen. I think eigen already does what you did there, no need to repeat it (?).
I mean you can pull a pointer, row/col count, + and * operator together and its a small matrix class, just call eigen for the * ? maybe that is useful to you, maybe not, hides the eigen part behind a simple interface... but again, I think their matrix already does that.

[baezacaljo]

Thanks jonnin, I got it. Also, I searched for some info, and now I have a better idea at the moment I wrote the previous code anything about pointers came to my mind, but after you mentioned, plus the small example in your previous reply. It seems I have a clue on thin wrapper, now I'll start implementing the information.

Thank you very much.