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 CONTINUEEND IFDO 20 i = 1, nrhs   CALL dtbsv( 'Upper', 'No transpose', 'Non-unit', n, kl+ku,$ ab, ldab, b( 1, i ), 1 )20 CONTINUEMy question is does a level 3 BLAS enabled, better performing variant of the forward/backward substitution routine exist?  
Level 3 BLAS accelerated solver for banded linear systems.
lapack;blas
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.