CODEWITHSUNDEEP
matrixlinear-algebralapack
enanone
enanone

Reputation: 966

LAPACK: Triangular system of equations with multiple right hand side

I have my beuatiful triangular n x n matrix, say L (for lower triangular), and I want to solve a system like

LX=B

Where B and X are n x k matrices (that is: I want to solve a triangular linear system with multiple right hand side). Additionally, I have my triangular matrix stored in PACKED FORMAT; i.e. I only store the lower triangular part. I am using BLAS and LAPACK, but I have realised that there is no specific way of solving my problem. Although there are many functions that solve similar problems:

stpsv(): Takes a triangular matrix in packed format and solves for a single right hand side.

strsm(): Takes a triangular matrix in dense format and solves for multiple right hand side.

What I really need is a combination of both. I would like to have a function accepting packed triangular format, as in stpsv(), and also accepting multiple right hand side, as in strsm(). But it seems as if there is not such a function readily available.

So my questions are:

  1. Is there any function that can accept a packed triangular matrix and solve for multiple right hand side?

  2. If the answer is NO, what would be more efficient? Either I call stpsv() in a for loop for every column in B, or I create a dense matrix from L, so that I have all those useless zeros in there and then I call strsm(). What would be better? Moreover, maybe I am missing a more clever way of doing all this.

Upvotes: 1

Views: 941

Answers (2)

francis
francis

Reputation: 9817

Yes, there is a function to solve Ax=B, for packed triangular matrix A and multiple right hand size B. It is stptrs() from LAPACK. In addition, there are other routines for triangular packed matrices, all featuring tp in their name according to the naming conventions of LAPACK.

However, looking at the source unveils that this function calls stpsv() from BLAS in a loop, once for each right hand side. It's exactly what you suggested!

Upvotes: 2

percusse
percusse

Reputation: 3106

Packed storage implies BLAS2 routines. Otherwise BLAS3 functions are more efficient in solving linear systems but they work in optimal blocked algorithms. If you call BLAS2 functions then you basically go back to vectored version hence it won't make too much sense.

Note that BLAS2 versions also do not perform conditioning checks. So they are directly optimized for BLAS2 performance since a triangular matrix with single RHS is a direct backward substitution.

For multiple RHS you can convert your matrix via, say stpttr and then use strtrs.

Upvotes: 3

Related Questions