Copying upper MatrixXd to lower MatrixXd (Eigen3) C++ library

Question

I've got a lower triangular MatrixXd and I want to copy its lower values to the upper side as it'll become a symmetric matrix. How can I do it?

So far I've done:

 MatrixXd m(n,n); 
 .....
 //do something with m
 for(j=0; j < n; j++)
       {
         for(i=0; i<j; i++)
           {
             m(i,j) = m(j,i);

           }
       }

Is there a fastest way to do it? I was thinking of some internal method that is able to "copy" the lower triangular matrix to the upper. Say I've got this matrix, we call m:

1 2 3
4 5 6
7 8 9

what I need to obtain in m is :

1 4 7
4 5 8
7 8 9

I also know you can get the upper or the lower part of the matrix to do something:

MatrixXd m1(n,n);
 m1 = m.triangularView<Eigen::Upper>();
cout << m1 <<endl;

1 2 3
0 5 6
0 0 9

But I can't yet get what I want...

Cusk

wht's the dimensions of your matrix? – Poleyn

@ Heisenbug The dimension is n (from 9000 to 32000) – Cusk

A little off-topic: is this C? Shouldn't it be m[j][i] = m[i][j]? – Muscular

@ EltanT m is MatrixXd object. The way to access to its members is defined that way. m is not a normal array – Cusk

@Cusk Okay. Anyway, you don't need to copy the diagonal, so the inner loop should run while i < j (and not i = n). Not a significant speedup I'm afraid, though. – Muscular

@ EltanT you are right there though – Cusk

Answer 1

I assume here that you are referring to working with the Eigen3 c++ library. This is not clear from your question. if not, you should consider it. In any case, within Eigen, there is no need to actually copy the triangular part, to get a selfadjoint matrix. Eigen has the concept of views, and you can use a self adjoint view in order to perform an operation like e.g.

using namespace Eigen;
MatrixXd m(m,n);
...
(generate uppper triangular entries in m)
...
VectorXd r(n), p(n);
r = m.selfadjointView<Upper>() * p;

here is a small example to illustrate using fixed size matrices:

#include <Eigen/Core>

using namespace std;
using namespace Eigen;

int main()
{
    Matrix2d m,c;
    m << 1, 2,
         0, 1;

    Vector2d x(0,2), r;

    // perform copy operation 
    c = m.selfadjointView<Upper>(); 
    cout << c << endl;

    // directly apply selfadjoint view in matrix operation
    // (no entries are copied)
    r = m.selfadjointView<Upper>() * x;
} 

the output will be [1, 2, 2, 1]. now, the result in r is the same as if you had used c * x instead. Just that there is no need for copying the values in the original matrix to make it selfadjoint.

Inhabited

Thanks for clarify. Yes, I meant the Eigen3 c++ library. However, I can't do what you say. First I've got class Eigen::MatrixXd' has no member named 'selfAdjointView from the compiler. On the other hand, what I need is to have the same elements on the upper and the lower side of the matrix. Therefore I need somehow to "copy" the elements from one side to the other within the same MatrixXd. – Cusk

I didn't compile the code, and there was a typo in selfadjointView. Fixed the code and added some clarifications. – Inhabited

Good! That works. I need to have all the values in the original matrix, m following your example. So what I've done is m=m.selfadjointView<Upper>();. However, it is slower than copying it with OpenMP. So I'll rather stick on my implementation. Thanks anyway! – Cusk

Answer 2

In case the selfadjointView is not an option for you, the solution is to use triangularView on the destination matrix:

m.triangularView<Lower>() = m.transpose();

Laryngoscope

Answer 3

The simplest way I can think of is by copying the upper part of m matrix trasposed on the upper part:

    m.triangularView<Upper>() = m.transpose();

For example, the following code:

    MatrixXd m(3,3);
    m << 1, 2, 3, 4, 5, 6, 7, 8, 9;

    m.triangularView<Upper>() = m.transpose();
    std::cout << m << std::endl;

Gives the output you asked for:

1 4 7
4 5 8
7 8 9

Regards.

Calliopsis

Answer 4

Simply:

m = m.selfadjointView<Upper>();

Leader

Answer 5

I think you are doing it the right way. If you knew some details about the memory layout of data in the matrix you could use some low-level optimizations. One of the techniques is loop tiling.

Hedgepeth

I am using OpenMP already, but I was thinking in some internal method of Eigen matrix that could speedup my code. – Cusk

I have seen a few presentations by Intel where they say they use tiling and OpenMP at the same time. – Hedgepeth

Answer 6

If speed is a big issue, I would not copy anything just decorate/wrap the matrix object with a coordinate inverting object that would flip the (x,y) to (y,x). if you make the () operator an an inline function it will not incur any significant cost when you use it.

Greatnephew

Answer 7

This works, you can cut something but you need at least n*m/2 (less something), so only of a 2x

edit: I see that you use this matrixd object... the syntax is different, but the algorithm is this, anyway

#include <stdio.h>


int main ( )
{
    int mat [ 4 ] [ 4 ];
    int i, j;

    mat [ 0 ] [ 0 ] = 0;
    mat [ 0 ] [ 1 ] = 1;
    mat [ 0 ] [ 2 ] = 2;
    mat [ 0 ] [ 3 ] = 3;
    mat [ 1 ] [ 0 ] = 4;
    mat [ 1 ] [ 1 ] = 5;
    mat [ 1 ] [ 2 ] = 6;
    mat [ 1 ] [ 3 ] = 7;
    mat [ 2 ] [ 0 ] = 8;
    mat [ 2 ] [ 1 ] = 9;
    mat [ 2 ] [ 2 ] = 10;
    mat [ 2 ] [ 3 ] = 11;
    mat [ 3 ] [ 0 ] = 12;
    mat [ 3 ] [ 1 ] = 13;
    mat [ 3 ] [ 2 ] = 14;
    mat [ 3 ] [ 3 ] = 15;

    for ( i = 0; i < 4; i++ )
    {
        for ( j = 0; j < 4; j++ )
            printf ( "%02d", mat [ i ] [ j ] );
        printf ( "\n" );
    }
    printf ( "\n" );

    for ( i = 1; i < 4; i++ )
    {
        for ( j = 0; j < i; j++ )
            mat [ j ] [ i ] = mat [ i ] [ j ];
    }

    for ( i = 0; i < 4; i++ )
    {
        for ( j = 0; j < 4; j++ )
            printf ( "%02d ", mat [ i ] [ j ] );
        printf ( "\n" );
    }
    printf ( "\n" );

    scanf ( "%d", &i );
}

Blithe

OP is using Eigen, not plain C. – Guyton

@Blithe Thanks for your help, but it is not that case. As I mentioned before, this is not a simple array in C. MatrixXd is an Eigen library object. – Cusk

Yes, but this does not change much: for(j=1; j < n; j++) { for(i=0; i<=j; i++) { m(i,j) = m(j,i); } } – Blithe

@mauropellizzer: That's exactly the OP's code, and OP wants something more efficient than this. – Guyton

No this not start from zero :) – Blithe