Level 3 BLAS accelerated solver for banded linear systems.

Question

At the moment I consider the following problem. I have a huge dense banded matrix $A$ which I want to factorize and use to solve linear systems $Ax=b$. $b$ has around more than 100 columns. At the moment I use DGBTRF from LAPACK to factorize $A$. This routine is a level-3 BLAS accelerated one which works fast and efficiently. The final step, solving $ Ax = LUx = b$ is done using DGBTRS which is not a level-3 BLAS accelerated one. In the current LAPACK version 3.5 this solver still works nearly sequentially on the columns of $b$. Namely the code looks like this:

IF( lnoti ) THEN
  DO 10 j = 1, n - 1
    lm = min( kl, n-j )
    l = ipiv( j )
    IF( l.NE.j ) CALL dswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
    CALL dger( lm, nrhs, -one, ab( kd+1, j ), 1, b( j, 1 ), ldb, b( j+1, 1 ), ldb )
10 CONTINUE
END IF
DO 20 i = 1, nrhs
   CALL dtbsv( 'Upper', 'No transpose', 'Non-unit', n, kl+ku,$ ab, ldab, b( 1, i ), 1 )
20 CONTINUE

My question is does a level 3 BLAS enabled, better performing variant of the forward/backward substitution routine exist?

Answer 1

Unfortunately, it looks like presently the only BLAS routines which take advantage of triangular/band structure are both level-2, _tbmv and _tbsv, the latter appearing in the code snippet above. (You are asking for a "_tbsm".)

The answer to your question depends on the interpretation of 'level-3 BLAS enabled'. There is presently no "_tbsm", so No. On the other hand, you could build "_tbmm" and "_tbsm" from the existing level-3 routines _trmm, _trsm, and _gemm; so Yes.

---

By the way, if this does turn out to be a missing feature, feel free to post a feature request in the LAPACK forums.