In case B(size ~ 2k, complex double) is block-diagonal, where block size is small(e.g. 2), is there any more efficient way to compute this other than Lapack gesv?

  • 1
    $\begingroup$ Assuming you already have an LU factorization, during the forward substitution step some elements of $B$ will be guaranteed zeros, so you can carefully invoke trsm to pass it only the relevant portion of the matrices for each block. This would save $\frac12$ of the first substitution step. $\endgroup$ – Kirill May 4 '17 at 17:13
  • $\begingroup$ @Kirill thanks, can you leave with an more explicit answer below(perhaps with psudo lapack flow)? suppose A is general and B is tridiagonal or 2x2block-diagonal. $\endgroup$ – lorniper May 7 '17 at 6:33

Supposing you already have an LU factorization, you can save a half of a forward substitution step. In the system $Lx=b$, you would have $$ \begin{pmatrix} l_{11}&0&&&&\cdots&0\\ \vdots&&&&&&0\\ l_{k-1,1}&\cdots&l_{k-1,k-1}&0&\cdots&\cdots&0\\ l_{k,1}&\cdots&l_{k,k-1}&l_{k,k}&0&\cdots&0\\ \vdots&&&&&&\vdots\\ l_{n1}&\cdots&&&&&l_{nn} \end{pmatrix} \begin{pmatrix} x_1\\\vdots\\x_{k-1}\\x_k\\\vdots\\x_n\end{pmatrix} = \begin{pmatrix} 0\\\vdots\\0\\b_k\\\vdots\\b_n \end{pmatrix} $$ and the observation is that only the submatrix $L[k:n,k:n]$ matters whenever the rhs vector $b$ has $b[1:k-1]=0$. (It makes no difference if only $b_k$ is nonzero, or if $b_{k:n}$ are nonzero.) The elements of $x$ corresponding to the first $k-1$ zero elements of $b$ all vanish.

So instead of invoking the usual procedure (in blas, this is trsm—you shouldn't have to write this yourself, and if you do it may well be slower) on the whole matrix, you invoke it only on the bottom-right submatrix of the right size, because you know that $x_1=\cdots=x_{k-1}=0$. For a $2\times2$ block-diagonal matrix, you can solve for two rhs $b$ vectors at a time with the same position of the first nonzero. This saves on average $\frac12$ of the work, averaged over all the $b$'s, and I believe this is the best you can do.

| cite | improve this answer | |
  • $\begingroup$ thanks so much! so does this also imply, in the case of diagonal right hand side, standard gesv is inefficient? $\endgroup$ – lorniper May 8 '17 at 13:52
  • $\begingroup$ One generally doesn't compute the matrix inverse explicitly, which is what one gets with a diagonal r.h.s., instead keeping the LU factorization itself, so it's basically never an issue. Life Stefano M points out, unless you've checked, and solving $A\backslash B$ is a major bottleneck, there's no real point in doing any of this in the first place. $\endgroup$ – Kirill May 8 '17 at 16:28

When solving $AX = B$ with dense $A$, the fact that $B$ is sparse does not imply that $X$ is sparse: this means that when computing $X$ little can be gained from the structure/sparsity of $B$.

It is true (as noted by Kirill) that with a custom implementation of the back solve step in the forward substitution step some floating point operations can be saved (since you know that some entries in the intermediate result are zero) but this implies introducing a for loop and multiple calls to trsm. I'm still convinced that a single call to trsm on the whole $L$ and $B$ compares favourably to multiple calls to trsm on sub-blocks of $L$ and $B$, unless you carefully optimise the block size of the custom implementation. (And of course, you have to take into account pivoting...)

Bottom line: you can save some operations, but to gain efficiency the implementation is still (in my opinion) non-trivial.

| cite | improve this answer | |
  • $\begingroup$ Back/forward substitution is a standard function in BLAS, you definitely don't really need to write your own: dtrsm. And you can invoke it so that it needs to work only on a portion of the matrices (using lda, ldb, the usual way). $\endgroup$ – Kirill May 7 '17 at 2:05
  • $\begingroup$ @Kirill I edited my answer for clarifying my point of view. Maybe my answer is not very informative, but I think that implementing your correct solution strategy is not so straightforward, on modern multi-core architectures. $\endgroup$ – Stefano M May 7 '17 at 22:51
  • $\begingroup$ Perhaps? I'm not sure how straightforward it is, my sense is that saving half the work is worth it, but, now that you mention it, I guess I wouldn't want to be too sure about it without measuring it first. (I haven't measured this.) $\endgroup$ – Kirill May 7 '17 at 23:25
  • $\begingroup$ @Kirill yes benchmarking is paramount: as this example shows, sometimes for loops can be even faster than equivalent vectorised code. Taking into account that you have LU factorisation, forward solve, and back solve, and that you save just half of the work on a single step, the maximum speedup you can expect is about 20% (or 6/5), if $A$ and $B$ have the same size. But again results depend on memory layout, processor architecture, compiler-level optimisations, so unless X=A\B is a major bottleneck I would not embark in your solution strategy. $\endgroup$ – Stefano M May 8 '17 at 7:54
  • $\begingroup$ I understand what you mean, but I don't know OP's situation, and there isn't enough information in the question itself to know if it's really worth doing, so I just answered the question as posted. Still, a saving of 20%, especially orthogonal to compiler optimizations and the like, is substantial enough that one might legitimately want to know about it: it doesn't strike me as totally unreasonable, like some microoptimization might. $\endgroup$ – Kirill May 8 '17 at 16:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.