3 ms·

Very true, and I think this is generally why C 'does it right'. But there are caveats. Firstly, almost no-one actually creates their matrices as literals. Inst

by rocqua 4y ago

Very true, and I think this is generally why C 'does it right'. But there are caveats. Firstly, almost no-one actually creates their matrices as literals. Instead they are build dynamically an usually incrementally. Nevertheless, people still think back to these very basic examples when trying to understand things, so this probably should be standard. Secondly, these decisions matter a lot for matrix vector multiplication. The 'C' way is very efficient for computing Ax which presumes vector x to be a Column vector. When you want to compute xA instead which presumes vector x to be a row vector, you want column-major ordering for efficiency. The same works for matrix-matrix multiplication AB where you want A to be row-major and B to be column-major. To sum up. Because of our writing direction, and the natural way to write 2-d arrays, it makes sense to have row-major ordering. That immediately means it also makes sense to use column vectors, generally speaking. There are also exceptions if you happen to know if a matrix will often be at the 'left' or 'right' side of a multiplication. Some more interesting notes that are less relevant: In C the matrix you defined has an incompatible type with how most matrices are stored (most people in C store rows contiguously, then define a matrix as an array of pointers to rows). int mat[3][3]; vs int *mat[3] These are meaningfully different. In the first case getting element `mat[1][1]` is just 'base-pointer + 4'. Whereas in the second case, getting element `mat[1][1]` is 'read pointer at base-pointer + 1 and add 1 to the added pointer`. Another common ways to use matrices is element-wise, where you refer to x[i][j]. Mathematicians expect i to be the row and j to be the column. Which matches C. Whilst Fortran does it the other way around.

bruce343434 4y ago
vs int (*mat)[3] = calloc(3, sizeof(int[3])); which is a pointer to (or array of) arrays of 3 ints. Laid out flat, and you still get the `mat[i][j]` notation.
- rocqua 4y ago
  That is indeed almost the same as int mat[3][3] The biggest reason people use the separate array of pointers is for very quick row swaps. Often comes up in matrix operations when trying to solve systems or put matrices in triangular forms.
owlbite 4y ago
> The 'C' way is very efficient for computing Ax which presumes vector x to be a Column vector. When you want to compute xA instead which presumes vector x to be a row vector, you want column-major ordering for efficiency. You have this the wrong way around. The C row-major ordering is better for A^Tx, the Fortran column-major format is better for Ax. If you do Ax in row-major, you end up multiplying a row of A by x as your simd vector op. This then leaves you needing to horizontally reduce the result, which you avoid by using a column-major layout that instead accumulates multiple different results.
- bruce343434 4y ago
  wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what? Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage?
  - kergonath 4y ago
    > wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what? No, the other way around. In row-major storage, all the elements of the same row are contiguous, so it is faster to process the array row by row. > Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage? The explanations here are decent enough. All you need is the figure at the top, really: https://en.m.wikipedia.org/wiki/Row-_and_column-major_order https://en.m.wikipedia.org/wiki/Row-_and_column-major_order .
  - adrian_b 4y ago
    The schoolbook formulae for most linear algebra operations are written using scalar products, e.g. for Ax as scalar products of rows of A with x. On a computer, such operations should almost never be implemented with scalar products, because the speed of a scalar product is limited by the latency of the FMA operation, not by its throughput. It is possible to accelerate the computation of a scalar product by ordering the elementary operations in a tree, but that complicates the program and it still does not reach the maximum speed of the hardware. Except for vector-vector operations (i.e. for level 2 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by AXPY operations (whose result vectors should be kept in registers, to not be limited by the memory transfer throughput), which are not limited by the latency of the FMA, like the scalar products. Except for vector-vector operations and matrix-vector operations (i.e. for level 3 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by rank-one matrix updates. Most of the computation of a rank-one matrix update consists of the tensor product of 2 vectors. Therefore, a matrix-matrix multiplication can be computed either by scalar products of the rows of the 1st matrix with the columns of the 2nd matrix, or by the tensor products of the columns of the 1st matrix with the rows of the 2nd matrix. The second method is much faster. The result of the tensor product must be kept in registers for the duration of the computation, so large matrices are partitioned in small blocks, i.e. submatrices, which are multiplied directly (an additional complication that increases the number of nested loops is that the small blocks that can be multiplied directly must be grouped in larger blocks that can be kept in the level-2 cache memory and reused for computations). For matrix-matrix multiplications it is less important whether they are stored in column major order or in row major order, because when submatrices are brought into the L2 cache, the matrix elements will be gathered or scattered to be in the best order, before being loaded into registers repeatedly. For matrix-vector multiplications, it is better to have the matrix in column major order, in order to compute the product by AXPY operations between columns and elements of the vector, and not by scalar products between rows and the vector (computing AXPY operations with shorter parts of a column at a time, so that the result vector can be kept in registers).
kergonath 4y ago
> Another common ways to use matrices is element-wise, where you refer to x[i][j]. Mathematicians expect i to be the row and j to be the column. Which matches C. Whilst Fortran does it the other way around. That’s not quite it. The fact that the row index is first is purely a convention. And by convention, i in In X(i,j) in some Fortran code is still the row index. The difference with C is the layout in memory, not the semantics: the fact that X(i+1,j) is contiguous to X(i,j) in Fortran and not in C, which means that for efficiency reasons inner loops in Fortran are over row indices rather than column indices as in C. It would have been odd for Fortran, which has been since the beginning a language for numerical calculations, to “not match what mathematicians expect”.

by rocqua 4y ago

bruce343434 4y ago

vs int (*mat)[3] = calloc(3, sizeof(int[3])); which is a pointer to (or array of) arrays of 3 ints. Laid out flat, and you still get the `mat[i][j]` notation.

rocqua 4y ago
That is indeed almost the same as int mat[3][3] The biggest reason people use the separate array of pointers is for very quick row swaps. Often comes up in matrix operations when trying to solve systems or put matrices in triangular forms.

owlbite 4y ago

> The 'C' way is very efficient for computing Ax which presumes vector x to be a Column vector. When you want to compute xA instead which presumes vector x to be a row vector, you want column-major ordering for efficiency. You have this the wrong way around. The C row-major ordering is better for A^Tx, the Fortran column-major format is better for Ax. If you do Ax in row-major, you end up multiplying a row of A by x as your simd vector op. This then leaves you needing to horizontally reduce the result, which you avoid by using a column-major layout that instead accumulates multiple different results.

bruce343434 4y ago
wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what? Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage?
- kergonath 4y ago
  > wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what? No, the other way around. In row-major storage, all the elements of the same row are contiguous, so it is faster to process the array row by row. > Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage? The explanations here are decent enough. All you need is the figure at the top, really: https://en.m.wikipedia.org/wiki/Row-_and_column-major_order https://en.m.wikipedia.org/wiki/Row-_and_column-major_order .
- adrian_b 4y ago
  The schoolbook formulae for most linear algebra operations are written using scalar products, e.g. for Ax as scalar products of rows of A with x. On a computer, such operations should almost never be implemented with scalar products, because the speed of a scalar product is limited by the latency of the FMA operation, not by its throughput. It is possible to accelerate the computation of a scalar product by ordering the elementary operations in a tree, but that complicates the program and it still does not reach the maximum speed of the hardware. Except for vector-vector operations (i.e. for level 2 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by AXPY operations (whose result vectors should be kept in registers, to not be limited by the memory transfer throughput), which are not limited by the latency of the FMA, like the scalar products. Except for vector-vector operations and matrix-vector operations (i.e. for level 3 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by rank-one matrix updates. Most of the computation of a rank-one matrix update consists of the tensor product of 2 vectors. Therefore, a matrix-matrix multiplication can be computed either by scalar products of the rows of the 1st matrix with the columns of the 2nd matrix, or by the tensor products of the columns of the 1st matrix with the rows of the 2nd matrix. The second method is much faster. The result of the tensor product must be kept in registers for the duration of the computation, so large matrices are partitioned in small blocks, i.e. submatrices, which are multiplied directly (an additional complication that increases the number of nested loops is that the small blocks that can be multiplied directly must be grouped in larger blocks that can be kept in the level-2 cache memory and reused for computations). For matrix-matrix multiplications it is less important whether they are stored in column major order or in row major order, because when submatrices are brought into the L2 cache, the matrix elements will be gathered or scattered to be in the best order, before being loaded into registers repeatedly. For matrix-vector multiplications, it is better to have the matrix in column major order, in order to compute the product by AXPY operations between columns and elements of the vector, and not by scalar products between rows and the vector (computing AXPY operations with shorter parts of a column at a time, so that the result vector can be kept in registers).

kergonath 4y ago

> Another common ways to use matrices is element-wise, where you refer to x[i][j]. Mathematicians expect i to be the row and j to be the column. Which matches C. Whilst Fortran does it the other way around. That’s not quite it. The fact that the row index is first is purely a convention. And by convention, i in In X(i,j) in some Fortran code is still the row index. The difference with C is the layout in memory, not the semantics: the fact that X(i+1,j) is contiguous to X(i,j) in Fortran and not in C, which means that for efficiency reasons inner loops in Fortran are over row indices rather than column indices as in C. It would have been odd for Fortran, which has been since the beginning a language for numerical calculations, to “not match what mathematicians expect”.