DMatrix multiplications produce wrong result

[seanxnie]

Hi, there I tried matrix multiplications with DMatrix and produced the wrong result. Here is a simple code for testing.

use nalgebra::{DMatrix}; // 0.33.0

fn main() {
    let nc: usize = 2;
    // upper triangle matrix filled by 1.0
    let a_mat = DMatrix::from_fn(nc, nc, |i, j| if i >= j {1.0 as f64} else {0.0 as f64}); 
    // diagonal matrix with diagonal elements filled starting from 1.0 with 1.0 increment.
    let p_mat = DMatrix::from_fn(nc, nc, |i, j| if i == j {i as f64 + 1.0} else {0.0 as f64}); 
    let t_mat = &a_mat * &p_mat;

    println!("A: {:?}\nP: {:?}\nT=A*P: {:?}\n", a_mat, p_mat, t_mat);
}

The output:

A: VecStorage { data: [1.0, 1.0, 0.0, 1.0], nrows: Dyn(2), ncols: Dyn(2) }
P: VecStorage { data: [1.0, 0.0, 0.0, 2.0], nrows: Dyn(2), ncols: Dyn(2) }
T=A*P: VecStorage { data: [1.0, 1.0, 0.0, 2.0], nrows: Dyn(2), ncols: Dyn(2) }

The result of T was supposed to be [1.0 2.0 0.0 2.0]. However, if we swap matrix A and P in this operation, ie, T = P*A, the result is correct. Therefore, I suspect the implementation of matrix multiplication is reversed. If this the case, kindly request to correct this issue. Thanks.

[sebcrozet]

Hi! The matrix multiplication result is correct. Note that nalgebra stores matrices in column-major format, so the result printed as a flat array might look transposed to you if you are used to row-major storage. I recommend printing the matrices with {} instead of {:?} to get a 2D view of the matrix array.

[seanxnie]

Thank you very much for the reply. As a newbie, I didn't notice that nalgebra's matrix is column-major based. If possible, please put more emphasis on the column-major underlying data in the documentation since it's a bit contradictory to commonsense. I appreciate that.

[Andlon]

@seanxnie: you're completely right, this is very confusing, and we need to address it. The way to resolve this is for our Debug implementation to use the same syntax as the matrix! family of macros, namely to use semi-colons to delimit rows. For example:

[ 1, 2, 3;
  4, 5, 6 ]

This way, it's unambiguous what the output means, and you can copy the output directly into a matrix! (or dmatrix!) invocation.

There's been work in #1119 to try to make Debug and Display output more consistent across many nalgebra types, but unfortunately the work stalled and nobody's been able to pick it up.