Discussion - Arrays vs Vector of Vectors when creating matrices - what is the most practical option?

Question

So i have stumbled upon an exercise on an C++ book called - "Guide To Scientific Computing in C++"

Here is the exercise: "Write code that dynamically allocates memory for three 2 × 2 matrices of double precision floating point numbers, A, B, C, and assigns values to the entries of A and B. Let C = A + B. Extend your code so that it calculates the entries of C, and then prints the entries of C to screen. Finally, de-allocate memory. Again, check you have de-allocated memory correctly by using a for loop as in the previous exercise."

This caught my eye, i tried to do this problem using both 2D Arrays (i was able to do it very easily) and with vector of vectors (i failed).

I did a lot of research and read some post on StackOverFlow and basically the opinion was universal - when dealing with matrices the option is always 2D arrays.

But knowing that there are a lot of programmers here (and that i am a newbie with C++), i would really like to read more opinions about this topic!

P.S: Here is a little snip of my failed attempt at trying to create a matrix using vector of vectors:

for (int row{ 0 }; row < 2; row++) { // Create Matrix A - goes through the matrix rows
    
    for (int col{ 0 }; col < 2; col++) { // goes through the matrix columns
        temp.push_back(rand() % 201); // add random number to the temporary vector
    }

    matrixA.push_back(temp);
}

// Outputing 
for (int row{ 0 }; row < matrixA.size(); row++) { // goes through the matrix rows

    for (int col{ 0 }; col < matrixA.at(row).size(); col++) { // goes through the vectors inside matrixA
        cout << matrixA.at(row).at(col) << "\t";
    }
    cout << endl;
}

This is the output: enter image description here

Answer 1

The root of the problem is that C++ brings no built-in support for multidimensional arrays, and

A matrix is not an array of arrays

So even though one can have in C/C++ an array of an array as e.g. in double M[3][3] this is not a good way to represent a matrix. Even if in the compile time fixed case this may still be ok, the vector< vector<double> > way is definitely not good.

Why:

• The subscript syntax M[i][j] is still reflecting the fact that this is an array of arrays and not a true multidimensional array (M[i,j] is legal in C/C++ but means something completely different). The M[i][j] syntax also does not really match the mathematics usage. To access the i,j element of a matrix mathematicians usually write

  A_i,j

  and not

  (A_i)_j

  (which would correspond to M[i][j])

• Experience in mathematics and scientific computing shows that most of the time both dimensions are usually equally important. There is not really an outer and inner dimension. Slicing is often done by both row and column. More complex subscript patterns (such as submatrix or striding) are also quite common.

• Frequently one needs a flattening operation, i.e. to reorganize the matrix elements in a linear vector.

• 0 x N or M x 0 matrices are not realizable. This applies to the static case double A[N][M] but also to the vector< vector<double> > where at 0 x N cannot be represented. Proper handling the empty matrices is extremely important. A matrix subscript may generate an empty matrix as a legitimate case. Matrix multiplication on the other hand will always need matching inner dimensions. So 0xN * NxM is a valid case and should trigger no error.

• The vector< vector<double> > wastes memory. Also it is slower than it needs to be (think of a Nx2 matrix with N large. this has a lot of overhead). Same would apply to double **M in pure C.

• vector< vector<double> > is not compatible with external linear algebra libraries (such as LAPACK or else)

• vector< vector<double> > could be non-rectangular, i.e. one needs to make sure that all inner vectors have the same length.

A matrix is an entity of its own

This is how it is done frequently in linear algrbra packages. There are differences but essetially are all very similar.

• It is most reasonable to treat a multidimensional array is an entity of its own. It is not an array of an arrays or so, even though it might seem to be a good idea at first sight.

• Storage of data is always linear

• The access of an matrix element is mapped into this linear storage. Actually it is quite simple. To access the i, j element, use A[i*lda + j] with lda the length of the leading dimension. Not that this implies a certain ordering of the element in memory. This case is commonly refered to as C order. The opposite, FORTRAN order, would just swap the meaning of i and j, or in case of higher dimensionality, reverse the order of the indices.

• To be really generic, a nd-array consists of

  struct ndarray {
    double *buffer;
    size_t  size[DIMS];
    size_t strides[DIMS]
  };
  

  This could implement a owning (i.e. has ownership over the buffer) or a non-owning matrix. The latter may be important to implement zero-overhead submatrix or column/row subscripting.

  Access to the i,j element of a matrix (DIM == 2) is then

  buffer[ i*strides[0] + j*strides[1] ]
  

  This readily generalizes to the higher dimensions.

  Also transposing of the matrix just requires reversal of the strides array. No copying is needed.

  A row or column subscript, a submatrix, or getting a strided submatrix could be implemented as a zero-overhead operation. One just needs to fill in the strides array properly and set buffer and size accordingly.

  strides has to be initialized propery on matrix creation, depending on whether C or FORTRAN layout is to be used.