Hello. I am a long time Fortran and C user, so pardon my inexperience with the C++ programming language. Over the years, I have developed many structural FEA programs for singular applications, and I realize it's time to evolve to a language which promotes generic programming. Clearly, Eigen is the front runner for such an application, and I would like to thank its authors for taking time to develop a robust package.

I have been pondering the following problem for the past few days: is there a way to populate a 3n x 3n matrix which references a sub matrix? For example, I need something to the effect of (which is written in not-so-correct syntax):

Code: Select all

Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> *M;
Eigen::Matrix3d subM;

//
 
M->block<3,3>(i,i) = &subM;

void f(){
   // operate on subM
}

//   after f(), M.block<3,3>(i,i) now points to the contents of subM


To achieve a high degree of efficiency, I want to avoid copying the contents of smaller matrices into once larger matrix. Can this be done? If so, how? I found this related link <viewtopic.php?f=74&t=82377>, but it does not seem to help me.

Hi,

is your goal to reduce memory usage or simply to avoid the copies?

If the latter, then the very little time you will save by skipping the copies will be completely lost at the first operation you do on it (because of cache coherence, vectorization, etc.).

The goal is to reduce memory usage. Is it correct to assume that redundant copying wastes memory and makes the program less efficient? I've always thought those two go hand-in-hand. (keep in mind I am dealing with large sparse matrices. Later on, once I understand the Eigen fundamentals, I will consider SparseMatrix).

Ok, if the destination matrix is sparse then what you are looking for would be a "block sparse matrix". We currently don't offer that. Nevertheless, since the matrix is sparse anyway, I strongly recommend you to directly deal with SparseMatrix using the devel branch.

I agree; sparse matrix is the way to go. With that being said, I still need to educate myself on the intricacies of Eigen, and C++ for that fact. Below, I provide a snippet of code demonstrate what I am trying to do.

Code: Select all

Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> M(n,n); // n = 6
M.setZero();
   
double *a = &M.block<3,3>(0,0).coeffRef(0,0);
Eigen::Map<Eigen::Matrix3d> m(a);

m(0,0)=1;
m(1,0)=2;
m(2,0)=3;
m(0,1)=4;
m(1,1)=5;
m(2,1)=6;
m(0,2)=7;
m(1,2)=8;
m(2,2)=9;

// output of m:
// 1 4 7
// 2 5 8
// 3 6 9

// output of M:
// 1 7 0 0 0 0
// 2 8 0 0 0 0
// 3 9 0 0 0 0
// 4 0 0 0 0 0
// 5 0 0 0 0 0
// 6 0 0 0 0 0



Ideally, what I want is:

Code: Select all

M = [1 4 7 0 0 0
     2 6 8 0 0 0
     3 6 9 0 0 0
     0 0 0 0 0 0
     0 0 0 0 0 0
     0 0 0 0 0 0]


Code: Select all

M.block<3,3>(0,0) << 1, 4, 7,
                                    2, 5, 8,
                                    3, 6, 9;


or:

Block<MatrixXd,3,3> m(M.block<3,3>(0,0));

then:

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

or:

m(0,0 = 1;
...


However, none of these methods will work with a SparseMatrix because you cannot reference a sub matrix of a sparse matrix (would be very inefficient).

Yikes... still having problems, this time with:

Code: Select all

Block<MatrixXd,3,3> m(M.block<3,3>(0,0));


My computer is throwing a fit. The compiler is saying:

Code: Select all

'Eigen::Block<XprType,BlockRows,BlockCols,InnerPanel,HasDirectAccess>' to 'const Eigen::Block<XprType,BlockRows,BlockCols>' &

I am using Eigen 3.0.3. Should I be using the Dev branch?

Nevermind... I figured this out. What I want is a block matrix. In other words,

Code: Select all

M = [a b
     c d]


Doing...

Code: Select all

Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> M(6,6);
M.setZero();

Eigen::Block<Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>,3,3> a(M,0,0);
Eigen::Block<Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>,3,3> b(M,3,0);
Eigen::Block<Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>,3,3> c(M,0,3);
Eigen::Block<Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>,3,3> d(M,3,3);

a << 1 , 2 , 3 ,
     4 , 5 , 6 ,
     7 , 8 , 9 ;

d << 11 , 22 , 33 ,
     44 , 55 , 66 ,
     77 , 88 , 99 ;


Yields:

Code: Select all

M = [1 2 3 0  0  0
     4 5 6 0  0  0
     7 8 9 0  0  0
     0 0 0 11 22 33
     0 0 0 44 55 66
     0 0 0 77 88 99]


This is what I was looking for all along. Thank you for your help and inspiration, ggael.