# Why does Eigen allocate a temporary to evaluate A.noalias() = B.transpose()*C in parallel?

I wrote a program which iteratively transforms data using matrix multiplications. To minimize the number of large memory allocations, I use two roughly equal-sized std::vector<double>'s data and temp_data to store the current data and the new data respectively. For the matrix operations, I wrap each buffer into an Eigen::Map<Eigen::MatrixXd>, called data_matrix and temp_data_matrix respectively. Then, I need to compute the following operation:

temp_data_matrix.noalias() = data_matrix.transpose()*factor;


where factor is simply an Eigen::MatrixXd which is much smaller than data_matrix or temp_data_matrix. Thanks to .noalias(), this line does not allocate any large temporaries normally. However, when multi-threading the program suddenly allocates a temporary (roughly) equal to the size of temp_data! In practice this can increase memory usage by around 50% so I'm curious about a few things:

• Why does this only happen when multi-threading? What algorithmic aspect requires a temporary buffer on top of the output buffer?
• Why does this not happen with the line temp_data_matrix.noalias() = factor*data_matrix? After all, my operation seems more cache-efficient: because we transpose data_matrix, we simply need to compute the scalar product of each pair of columns from data_matrix and factor, which should be fast considering both matrices are column-major. However, I'm not considering the effects of blocking here, which I know can be very important for performance.
• What would be the fastest way to implement the aforementioned line without allocating any temporaries of the order of data or temp_data? In case the parallel algorithm needs some "scratch space", can't it use temp_data_matrix for that? I don't care about copying factor because it's relatively small in practice. Right now the fastest solution seems to be to switch back to single-threading for this line, but I hope this can be improved. I tried using lazyProduct but this turned out to be much slower than the single-threaded solution.

Appendix: The observations above were made using Massif (valgrind) both in my original program and with the example code below. I compiled the example program using the command g++ -O3 -DNDEBUG -I /usr/include/eigen3 -fopenmp eigen_memory_test.cpp -o eigen_memory_test.

#include <iostream>
#include <Eigen/Dense>
#include <algorithm>
#include <random>
#include <cstring>

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

bool parallel = (argc > 1 && strcmp(argv[1], "parallel") == 0);
std::cout << "parallel: " << (parallel  ? "true" : "false") << std::endl;
bool lazyproduct = (argc > 2 && strcmp(argv[2], "lazyproduct") == 0);
std::cout << "lazyproduct: " << (lazyproduct  ? "true" : "false") << std::endl;
bool test = true;

// Initialize parallellism
int threads = 4;
Eigen::setNbThreads(threads);

// Initialization and memory allocation
size_t columns = 100000;
std::vector<size_t> rows_per_step = {100, 90, 80, 70, 60, 50};
size_t steps = rows_per_step.size() - 1;
std::vector<double> data(rows_per_step[0]*columns);
std::vector<double> temp_data(rows_per_step[0]*columns);

// Generate random data
std::random_device rnd_device;
std::mt19937 mersenne_engine {rnd_device()};
std::uniform_real_distribution<double> distribution(-1.0, 1.0);
auto generator = std::bind(distribution, std::ref(mersenne_engine));
std::generate(data.begin(), data.end(), generator);

// If testing, print first elements to verify later
if (test)
std::cout << "data[:3]: [" << data[0] << ", " << data[1] << ", " << data[2] << "]" << std::endl;

// Process steps
for (size_t step = 0; step < steps; step++) {

size_t initial_rows = rows_per_step[step];
size_t new_rows = rows_per_step[step + 1];

// Choose factor
Eigen::MatrixXd factor = (test ?
Eigen::MatrixXd(Eigen::MatrixXd::Identity(initial_rows, new_rows)) :
Eigen::MatrixXd(Eigen::MatrixXd::Random(initial_rows, new_rows))
);

// Transform data
Eigen::Map<Eigen::MatrixXd> data_matrix(data.data(), initial_rows, columns);
Eigen::Map<Eigen::MatrixXd> temp_data_matrix(temp_data.data(), columns, new_rows);
if (!parallel)
Eigen::setNbThreads(1);
// In my real problem setting, the transpose in the following matrix multiplication is actually necessary
if (lazyproduct)
temp_data_matrix.noalias() = data_matrix.transpose().lazyProduct(factor);
else
temp_data_matrix.noalias() = data_matrix.transpose()*factor;
Eigen::setNbThreads(threads);

// Transpose back to original data
// This doesn't happen in the real problem, but I put this here to keep this example simpler
Eigen::Map<Eigen::MatrixXd> transposed_data_matrix(data.data(), new_rows, columns);
transposed_data_matrix.noalias() = temp_data_matrix.transpose();

}

// Print first elements to make sure matrix operations aren't optimized away
std::cout << "data[:3]: [" << data[0] << ", " << data[1] << ", " << data[2] << "]" << std::endl;

return 0;

}
$$$$
`