parsec
parsec copied to clipboard
parsec_matrix_define_triangle swaps Upper/Lower Triangular matrices?
Original report by Omri Mor (Bitbucket: omrimor2, GitHub: omor1).
diag = /* 0 if includes diagonal, 1 otherwise */
if (uplo == matrix_Upper) {
nmax = n-diag;
for (i = diag; i < n; i++) {
unsigned int mm = i + 1 - diag;
blocklens[i] = mm < m ? mm : m;
indices[i] = i * ld;
}
} else if (uplo == matrix_Lower) {
nmax = n >= (m-diag) ? m-diag : n;
for (i = 0; i < nmax; i++) {
blocklens[i] = m - i - diag;
indices[i] = i * ld + i + diag;
}
diag = 0;
}
parsec_type_create_indexed(nmax, blocklens+diag, indices+diag, oldtype, &tmp);
Unless I’m misreading this or fundamentally misunderstanding how PaRSEC treats upper and lower triangular matrices, the definitions of the two here are swapped. An upper triangular matrix should have min(n, m)
rows, with the first row having m
elements, the next having m-1
, and so forth, until the last row has 1
element. A lower triangular matrix has n
rows, with the first have 1
element, the following 2
, and so forth, until it caps at min(n, m)
elements.
Unless PaRSEC uses column-major storage order, in which case n
is the number of columns, rather than rows, as I’d understood it?