I have an upper Hessenberg matrix i.e. upper trapezoidal with a block=1 band under the diagonal with non zero elements from some column onwards, and I need to triangularize it in the fastest possible way. This use-case is the result of computing e.g. QR and then deleting one column (yes it is always only one column) from R. I need to make R upper triangular again by annihilating all those (one per column) non-zero elements below the diagonal e.g. I need to get rid of the 2s in the fastest way:
- Code: Select all
1 1 1 1 1
0 1 1 1 1
0 2 1 1 1
0 0 2 1 1
0 0 0 2 1
I have tried:
- Using Givens rotations combining dlartg and dlasr in a loop. Issue: very bad spatial locality, keeps moving along the rows that consequences in high cache misses.
- Using householder reflectors combining dlarfg and dlarf in a double loop. Issue: very bad spatial locality, keeps moving along the rows that consequences in high cache misses.
- Using the latest 3.4.x QR decompositions e.g. dgeqrt. This is a work overkill for this simple update, algorithmically very bad even when the column deleted is among the firsts ones.
Can anyone point me to a set of lower level functions I can use to implement my use-case in the fastest way? Otherwise I might end up reinventing the wheel, studying the update pattern and doing it manually in a block of registers along the columns.