# Vectorize a part of a huge matrix in C++

I have a large matrix, side-length is about $n\geq 1000$. I need to do element-wise multiplication of this matrix with another matrix many, many times. I make this process by:

1. Vectorizing (through armadillo) this matrix, i.e., flattening it in a single array.
2. Splitting this array over many cores. Each core will take a part of this matrix and do its multiplication.

Now I realized that for Mathematical reasons, I only need the upper-triangular part of this matrix. How can I vectorize only the upper-triangular part and discard the other components of the matrix? Simply looping is not that good. Do I have any other options?

My matrices are Armadillo complex matrices of arma::cx_mat. And my current vectorization function is:

typedef double Real;
typedef std::vector<std::complex<Real> > VecType;
arma::cx_mat myMatrix;
//fill my matrix...
VecType  vectorizedMatrix = arma::conv_to<VecType>::from(myMatrix);


I find std::vector very useful, because I do element-wise multiplication very efficiently using standard algorithms with std::transform(), which is best optimized (AFAIK).

• Are you able to avoid creating the lower triangular part in the first place? I'm sure there are some matrix libraries that store just the upper-triangular part of a matrix as a single array, and then "vectorizing" the matrix is a no-op (you'll have to be sure to check row-major vs column-major storage). It's hard to imagine doing better than that. (Also, I believe most "symmetric" matrix types in matrix libraries are only storing a triangle of the matrix, so you might be able to trick that into doing what you need, even if your problem isn't symmetric.) – Tyler Olsen Jul 18 '16 at 16:35
• @TylerOlsen Well, you're right... but I use Armadillo... It would be more expensive to convert all matrices back and forth. – The Quantum Physicist Jul 18 '16 at 17:00
• Why don't you want to just loop through the upper triangular part of the matrix once and manually fill in vectorizedMatrix? If you first reserve() or resize() with the right amount of space this should be very fast. Also be sure to properly order your loops for cache efficiency. – Doug Lipinski Jul 18 '16 at 22:52
• @DougLipinsksi That's my last resort actually. I checked Armadillo's source code internally and found that it does something similar but with some memory copying, and not some simple assignment operator. That's probably why it's fast. I'll try to do that. – The Quantum Physicist Jul 19 '16 at 7:34
• Just to have one more google search term: the operation that you need is known as half-vectorization. – Federico Poloni Aug 21 '16 at 10:29

So I invented my own solution. Critique, suggestions are highly appreciated.

I learned that elements are stored in a matrix as column-major. So I assumed that using std::copy() with an iterator over a single column will just give consecutive elements, hence better speed.

Here's my solution:

template <typename T>
std::vector<T> vectorizeMatrix_upperTriangular_columnMajor(const arma::Mat<T>& mat)
{
std::vector<T> output((mat.n_cols*(mat.n_cols-1))/2);
typename arma::Mat<T>::const_iterator it = mat.begin() + mat.n_rows; //iterator at rows to skip in every step, starts at second column
long toSkipInVec = 0;
for(int i = 1; i < mat.n_cols; i++) //Starts with 1 to skip the diagonal
{
std::copy(it, it + i, output.begin() + toSkipInVec);
toSkipInVec += i;
it += mat.n_rows;
}
return output;
}