Solving a linear system with triangular matrix

[Mikhail Katliar]

I want to solve the system of equations A*x=b, where A is lower-triangular. I can use either the blas::trsv() function or the inv() expression:

blaze::LowerMatrix<blaze::DynamicMatrix<Real, blaze::columnMajor>> A(M, M);        
blaze::DynamicVector<Real, blaze::columnVector> b(M);
blaze::DynamicVector<Real, blaze::columnVector> x(M);

randomize(A);
randomize(b);

bool const use_inv = true;
if (use_inv)
  x = inv(A) * b;
else
{
  x = b;
  trsv(A, x, 'L', 'N', 'N');
}

• I noticed that for big A matrices the inv(A) * bversion results in a call to BLAS dtrtrifunction, i.e. the matrix is first inverted and then multiplied with a vector. This is less efficient than using trsv(). I think it would be an improvement if Blaze could detect that A is a LowerMatrix and use trsv()instead.
• I see from the code that for small ABlaze uses custom inversion kernels. But still, it first inverts A and then multiplies it with b.
• I could directly call trsv(), but I want to take the advantage of custom kernels for small matrices, which are, according to my benchmarks, significantly faster than calls to BLAS.

Would it be possible to make it such that when evaluating inv(A) * bwith triangular A

1. for large matrices, BLAStrsv()is used,
2. for small matrices, solution is performed using optimized small kernels,
3. no inversion of A is performed?

Best regards,

Misha

[Klaus Iglberger]

Hi Misha!

You are absolutely correct, in this case there should be no inversion, but only the back substitution. If there is indeed an inversion happening then this is an oversight/defect on our side. We will fix this as soon as possible. Thanks a lot for pointing this out,

Best regards,

Klaus!

[Mikhail Katliar]

Hello Klaus,

Is there an explicit way to tell Blaze to perform the substitution and not matrix inversion (besides using BLAS trsv wrapper)?

[Klaus Iglberger]

Hi Misha!

No, so far there is no way to explicitly perform a substitution (for instance in form of a solve() function). Your best option is to call trsv() explicitly.

Best regards,

Klaus!

[Klaus Iglberger]

The feature has been implemented, tested, optimized, and documented as required. The following four expressions are now handled by means of linear system solvers instead of inverting a matrix and a subsequent multiplication:

• x = inv(A) * b; where x and b are dense column vectors and A is a dense matrix;
• x = b * inv(A); where x and b are dense row vectors and A is a dense matrix;
• X = inv(A) * B; where X, B, and A are dense matrices;
• X = B * inv(A); where X, B, and A are dense matrices.

The implementation adheres to the three given requirements:

• for large matrices, BLAS trsv() or faster approaches are used;
• for small matrices, the solution is computed by using manually optimized kernels;
• no inversion of A is performed.

It is immediately available via cloning the Blaze repository and will be officially released in Blaze 3.7.

[Mikhail Katliar]

Great, thanks for the update!