Directly accessing entries of armadillo-matrix by memory-pointer (memptr) does not work after matrix is modified

I have a problem where I want to access certain entries of an armadillo-matrix "M" by a pointer in a struct (or class). After initializing M I set the pointer in the struct. By dereferencing the pointer I can see it has the right value (the first entry of M - or M(0,0)).

Then I change M to M * M. But now dereferencing the pointer does not give me the right value anymore. What's weird: If I have a small matrix i.e. 3x3 or 4x4 (see "matrixSize" in my code) the error does not happen. With small matrices dereferencing the pointer gives the RIGHT value. Bigger matrices though result in the wrong values (with 5x5 something like "2.76282e-320", which is probably some random place in memory).

What am I doing wrong here? How can I solve this problem?

If it helps, I'd like to explain what I want to achieve: I have a network of delay-coupled nodes that each have some sort of dynamic behaviour. (think delay-coupled differential equations DDEs - delay-coupled Oscillators). As they are delay-coupled I need to store their past states (an array of their histories). Each oscillator also has some LOCAL dynamic with dynamical variables that are not influencing other nodes which means that I don't have to store their histories. The matrix shall be used to keep the past states of some variable of the nodes. I want to have them in a matrix, because I want to use vector-operations on them (one index of the matrix represents time, while the other is the node-index). But I also want to access them individually to calculate some local dynamic at each of the nodes (oscillators). So I want to update the individual nodes, but also the global state. That's why having both representations helps: For the local dynamics I access the delayed states through a pointer into the matrix. For the global dynamics I access the coupled variables through a row in the matrix that functions as a history-array.

#include <iostream>
#include <armadillo>

struct someStruct {
    const double *currentStateInMatrix;
    void printvalue() {std::cout << "*currentState (pointer to first enty of matrix) = " << *currentStateInMatrix << std::endl;}
};

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

    uint M_size = 4;

    arma::Mat<double> M(M_size, M_size, arma::fill::randu);


    double* pointerIntoMatrix = M.memptr();
    std::cout << "pointer [ " << pointerIntoMatrix << " ] has value: " << *pointerIntoMatrix << std::endl;

    someStruct myStruct;
    myStruct.currentStateInMatrix = pointerIntoMatrix;

    std::cout << "original matrix M: \n" << M;
    myStruct.printvalue();
    std::cout << "\n+++changing contents of matrix M+++\n\n";

    M = M * M;

    std::cout << "M * M: \n" << M;
    myStruct.printvalue();
}

I compiled it using: g++ file.cpp -std=c++17 -larmadillo -o main

Solution

The operation M = M * M in Armadillo is a matrix multiply (not an element by element multiply). So storing the intermediate calculations of M * M directly into M would be problematic. It would overwrite existing data in M that is still needed to complete the M * M operation.

It's probably safe to assume that Armadillo detects this problem and stores the result of M * M into a separate chunk of memory and then assigns that chunk to M.

There are ways around that. Use fixed size matrices like darcamo mentioned in the comments. Use Mat<double>::fixed<4, 4> M; to declare the matrix.

Another way is to manually manage the memory for the matrix elements and tell the M matrix to always use that memory. There are the advanced constructors to do that:

mat(ptr_aux_mem, n_rows, n_cols, copy_aux_mem = true, strict = false)

So we can do:

double* my_mem = (double*)std::malloc(M_size * M_size * sizeof(double));
// or use new[]

Mat<double> M(my_mem, M_size, M_size, false, true);

// don't forget to call std::free(my_mem) or delete[] when done!

A word of caution. The above two workarounds (fixed size matrices and manual memory management) probably have minor performance penalties. There is no way to avoid using a separate chunk of memory for the M = M * M operation. When using these workarounds, Armadillo will store the result of M * M into a temporary buffer and then copy the contents of the buffer into the memory used by M.