blaze::trsv() returning incorrect result for row-major matrix argument

[Mikhail Katliar]

Example:

#include <blaze/Math.h>
#include <limits>
#include <cassert>

void test()
{
    // Change rowMajor to columnMajor to make it pass
    blaze::DynamicMatrix<double, blaze::rowMajor> A {
        {1., std::numeric_limits<double>::signaling_NaN()},
        {-2., 3.}
    };

    blaze::DynamicVector<double> x {4., 7.};
    trsv(A, x, 'L', 'N', 'N');
    blaze::DynamicVector<double> expected {4., 5.};

    assert(x == expected);   // <-- fails
}

[Mikhail Katliar]

This is the piece of code in blaze/math/lapack/trsv.h which is supposed to handle the row-major matrix case:

   if( IsRowMajorMatrix_v<MT> ) {
      ( uplo == 'L' )?( uplo = 'U' ):( uplo = 'L' );
   }

Passing a row-major matrix A to LAPACK's trmv() function which expects a column-major matrix effectively transposes it, changing a lower-triangular matrix to an upper-triangular one and vice versa. The resulting system being solved is then A^T x = b, which is different from the original one A x = b, hence the incorrect result.

For real matrices, this can be fixed by inverting the transpose flag for row-major matrices:

   if( IsRowMajorMatrix_v<MT> ) {
      ( uplo == 'L' )?( uplo = 'U' ):( uplo = 'L' );
      ( trans == 'N' )?( trans = 'T' ):( trans = 'N' );
   }

However, for complex matrices when trans='C' this will not produce the correct result, because changing trans from 'C' to 'N' will fix the element order but will not conjugate the elements as required. The system actually solved will be A^T x = b instead of A^H x = b.

[Klaus Iglberger]

Hi Mikhail!

I believe I understand the problem: You might have missed the specification for the LAPACK wrapper functions at the beginning of the Linear System Solver wiki entry:

[...] Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems:

Single right-hand side:

• A *x=b if A is column-major
• A^T*x=b if A is row-major

I hope that with this knowledge you can use the trsv() functions without problem. I apologize for the confusion, but there unfortunately seems no easier, but still consistent way to wrap the LAPACK solver functions.

Thanks for creating the issue (I admit this is a complex topic),

Best regards,

Klaus!