Error in evaluation of A*inv(trans(B)) expression

[Mikhail Katliar]

Code:

size_t constexpr M = 7;

StaticMatrix<double, M, M, columnMajor> L, B, X;
randomize(L);
randomize(B);

auto X1 = evaluate(B * inv(trans(L)));
auto X2 = evaluate(B * evaluate(inv(evaluate(trans(L)))));

std::cout << "X1 = \n" << X1 << std::endl;
std::cout << "X2 = \n" << X2 << std::endl;

std::cout << "X1 * trans(L) - B = \n" << X1 * trans(L) - B << std::endl;
std::cout << "X2 * trans(L) - B = \n" << X2 * trans(L) - B << std::endl;

Example output:

X1 = 
(     -1.28311     0.811128    -0.930073    -0.481951      1.26091      0.17977    -0.541146 )
(     0.165468     0.630649    0.0292309    -0.463089     0.549342     0.318625    0.0401133 )
(      0.16363     0.203228     0.582571     0.775679   -0.0829435      0.52525      1.00258 )
(     0.239505    -0.108393     0.938073     0.118108     0.704706    -0.265765   -0.0770531 )
(     0.923148    -0.182332     0.263774     0.891805    -0.615401     0.595352      1.92369 )
(      1.46814    -0.660995      0.65819      0.78423     -1.23826    -0.215314     0.175482 )
(    -0.237447     0.157807   -0.0205074     -0.54736     0.256302    -0.355022     -1.00695 )

X2 = 
(     -1.12063     0.552327     0.293299    -0.538178      1.39998      1.81349    -0.510475 )
(     -2.16823      0.76024     0.721346   -0.0677476      2.08721      1.56163    -0.247799 )
(     0.327226   -0.0263235     0.288595   -0.0283733    -0.259541     0.180702     0.253264 )
(    -0.773121     0.854442     0.470443    -0.187619      1.32012     0.425175    -0.624902 )
(    -0.488366     0.162849     0.123909     0.234593      1.20025      1.07009     -0.97149 )
(     -1.77946     0.942372     0.458295     0.487555      1.15168     0.618759    0.0689882 )
(     0.166745    -0.170762     0.306025    -0.453909     0.753515     0.837677    -0.734451 )

X1 * trans(L) - B = 
(     -2.30086     -1.20839       -2.341     -2.58472     -1.29465     -0.74623     -2.23135 )
(    0.0678443    -0.280896   -0.0429704      0.50211    -0.306903       0.1749    -0.547458 )
(      1.21319     0.809892      1.12819      1.13059      1.24864      1.27807      1.59177 )
(     0.821675    0.0295633     0.567464     0.486073    -0.395484     0.429763     0.769448 )
(      1.55505     0.410157      1.56235       2.6056      1.00916      2.31739      3.11656 )
(      2.03627   -0.0348916      1.38906      1.44798    -0.185528    -0.236153     0.568919 )
(     -1.63748    -0.580973     -1.98577     -1.72269     -1.07666     -1.27013     -1.73343 )

X2 * trans(L) - B = 
( -9.99201e-16 -4.44089e-16 -8.88178e-16 -3.33067e-16 -4.44089e-16 -1.33227e-15 -6.66134e-16 )
( -5.55112e-16 -3.33067e-16 -3.33067e-16 -5.55112e-17 -2.22045e-16 -9.99201e-16 -2.22045e-16 )
( -1.11022e-16            0            0            0 -5.55112e-17 -3.33067e-16  2.22045e-16 )
( -3.33067e-16 -2.22045e-16  5.55112e-17            0            0 -2.22045e-16 -2.77556e-16 )
( -5.55112e-16 -2.22045e-16 -5.55112e-16 -1.11022e-16 -2.22045e-16 -9.71445e-16 -2.35922e-16 )
( -6.03684e-16 -2.22045e-16 -2.11636e-16 -3.95517e-16 -1.11022e-16 -6.66134e-16 -1.11022e-16 )
( -2.22045e-16            0 -2.22045e-16 -1.11022e-16 -1.66533e-16 -3.60822e-16 -2.77556e-16 )

The X1 value is wrong, while X2 is correct. Interestingly, I could not reproduce this bug for M < 7.

[Klaus Iglberger]

Hi Misha!

Thanks a lot for pointing out this defect. You are correct, for N > 7 the combination of inv() and trans() unfortunately computes the solution for the transpose matrix. We apologize for the inconvenience and will fix the problem as quickly as possible.

Best regards,

Klaus!

[Klaus Iglberger]

Commit bdc64ea fixes the computation resulting from the expression A*inv(trans(B)) for column-major right-hand side matrices and system matrices of size N > 7. The fix is immediately available via cloning the Blaze repository and will be officially released in Blaze 3.7.