Updating Computing the T factor for OR/UNGQR family, - #1373
Conversation
…currently debugging
…currently debugging
…rft_ut version as a bailing of clarft
…the real and complex versions in how things are computed due to me not having the kernel to compute y=conjg(A)x+y for triangular y
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1373 +/- ##
==========================================
- Coverage 69.01% 68.66% -0.36%
==========================================
Files 6122 6166 +44
Lines 486123 489976 +3853
Branches 23286 23286
==========================================
+ Hits 335514 336418 +904
- Misses 150420 153369 +2949
Partials 189 189
... and 5 files with indirect coverage changes Continue to review full report in Codecov by Harness.
|
|
@jprhyne fyi, the failing extended-api-only tests, mean that some routines are missing external declarations for functions or subroutines that are used. For example, SSYTRRK uses SSYTRRK_LVL2 but never declares it. |
|
Thank you! I forgot to add those! I'll get that fixed and pushed shortly |
|
It is a beautiful pull request @jprhyne. Good job! I love it. ❤️ |
…on in the complex out of place trmm
|
@jprhyne I think you could use - CALL DSYTRRK('Lower', 'Upper', 'Transpose', 'Unit', K, ONE,
- $ V, LDV, ZERO, T, LDT)
+ DO J = 1, K
+ T(J,J) = ONE
+ DO I = 1, J-1
+ T(I,J) = V(J,I)
+ END DO
+ END DO
+ CALL DLAUUM('Upper', K, T, LDT, INFOThat way we could get rid of these new BLAS-like routines: As far as I can tell, we don't need those anywhere else. So adding them just to replicate this existing functionality seems a bit overkill and redundant to me. Performance-wise, it would be the same, since |
|
@jprhyne You can ignore the errors in the AppVeyor pipelines. They are caused by the recent LLVMFlang 23.1.0 release, which seems to break some Fortran standard module paths... |
|
@ACSimon33 Thanks for that catch! I planned to use the While these paths are not hit yet by any code (they will be in a future PR I want to make), this behavior is already present in If it is desired, I am more than happy to remove this behavior from just |
|
I feel that I should explain more about why I think For the
where In order to use However, the only upper triangular matrix available is This means that we would need to compute
|
|
@ACSimon33 I just noticed after I sent that previous comment, we also need this new routine for the In the
where Since we need to store the lower triangular component, we must enter
unless there's another routine that I've missed that does compute these two operations |
| * Base case | ||
| * | ||
| IF(N.EQ.1.OR.K.EQ.1) THEN | ||
| IF(K.EQ.1) THEN |
There was a problem hiding this comment.
Only this precision on the routine removed the N.EQ.1.OR. from the quick return (also see line 207: N.EQ.0.OR.). That can result in garbage output or errors:
DLARFT N=1, K=2 -> ** On entry to DSYRK parameter number 4 had an illegal value
DLARFT N=0, K=2 -> ** On entry to DSYRK parameter number 4 had an illegal value
| INFO = 11 | ||
| END IF | ||
| IF (INFO.NE.0) THEN | ||
| CALL XERBLA('DTRSM ',INFO) |
There was a problem hiding this comment.
XERBLA('DTRSM ') called from CTRSM_MOD. This is one of 15 sites where the name string does not match the routine:
{s,d,c,z}trsm_mod.f:234 XERBLA('DTRSM ')
{s,d,z}trsm_lvl2_mod.f:243 XERBLA('DTRSM ')
ctrsm_lvl2_mod.f:243 XERBLA('CTRSM ')
{s,d,c,z}trti2_mod.f:161 XERBLA('DTRTI2')
{s,c,z}trtri_mod.f:158 XERBLA('DTRTRI_MOD')
| * Definition: | ||
| * =========== | ||
| * | ||
| * SUBROUTINE DLARFT_UT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) |
There was a problem hiding this comment.
APPLYT missing. See other precisions as well.
| *> | ||
| *> \verbatim | ||
| *> | ||
| *> DLARFT forms the triangular factor T of a real block reflector H |
There was a problem hiding this comment.
The new DIRECT=T and STOREV=T options are not documented in {s,d,c,z}larft.f and {s,d,c,z}larft_ut.f. It would be good for users to know, since these are public API changes.
I built and benchmarked it, and I think the
|
| n | K | LVL2 |
PR NX=0 |
PR NX=64 |
xLAUUM |
xGEMV |
LVL2 (MKL) |
PR NX=0 (MKL) |
PR NX=64 (MKL) |
xLAUUM (MKL) |
xGEMV (MKL) |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 512 | 16 | 31.6 | 69.1 | 45.9 | 45.8 | 45.3 | 7.5 | 89.7 | 40.6 | 39.6 | 40.0 |
| 512 | 32 | 127 | 265 | 174 | 173 | 172 | 30.0 | 323 | 144 | 144 | 142 |
| 512 | 64 | 491 | 1024 | 664 | 655 | 653 | 116 | 1212 | 547 | 530 | 501 |
| 1024 | 16 | 67.0 | 109 | 86.0 | 85.9 | 85.4 | 15.3 | 99.0 | 50.0 | 50.2 | 48.3 |
| 1024 | 32 | 273 | 422 | 330 | 328 | 327 | 57.7 | 360 | 166 | 164 | 168 |
| 1024 | 64 | 1084 | 1634 | 1276 | 1266 | 1264 | 219 | 1286 | 590 | 581 | 576 |
| 2048 | 16 | 138 | 189 | 166 | 166 | 166 | 27.6 | 117 | 66.6 | 65.3 | 65.0 |
| 2048 | 32 | 563 | 731 | 640 | 638 | 637 | 108 | 398 | 214 | 214 | 211 |
| 2048 | 64 | 2266 | 2856 | 2498 | 2487 | 2486 | 433 | 1368 | 729 | 678 | 685 |
Once NX is set the three formulations are basically indistinguishable — 0.3–1.8% spread on minimums with the reference BLAS. Note the whole grid above is K <= 64, since ILAENV(3,'DLARFT') caps the UT path at 64 and larger blocks never reach it. Note also that ILAENV(1,'xLAUUM') is 64, so on a k-by-k block xLAUUM itself falls through to the unblocked xLAUU2. The blocked path is dead there too, which is why the xLAUUM and xGEMV variants tie (another reason to prefer the xGEMV version).
|
I don't really follow the block. If there was a bug that was traced back to this, I would really hate to be the one to debug this. In addition, it is a bit involved for what I would expect to be a reference implementation. If this is the I do like the second one, it is very clean! However for a later PR where we are computing a panel factorization, we will need to compute the following which will mean to compute which will require the routine So since there is two places that I have where we need One thing to note is also that your table shows the 3 methods are basically the same in execution time (which
One more note on point number 2, some experiments I have ran show that trmm is about twice as fast as trmmoop, indicating I am also interested in how this behaves when comparing against AOCL since that's been my main |
The reference BLAS are not useful for performance comparisons. Use any optimized BLAS -- OpenBLAS, MKL, BLIS, ... |
4 participants
Description
In order to make my ORGQR optimization PR more digestible and revertible, I am breaking up my original draft PR into different parts. This part is concerned with the computation of the triangular factor and adding the UT-based transform discussed
by Joffrain et. al: https://doi.org/10.1145/1141885.1141886
Following their notation, we will refer to$S$ as the triangular factor that must be applied with a solve and $T$ to be the triangular factor that must be applied with a multiplication (Note that $T=S^{-1}$ )
We added a few new blas-like kernels (listed first) that are not currently optimized by vendors but rely on existing routines where possible. These kernels are required for either using the new factor$S$ or computing the existing factor $T$ .
For all recursive kernels, we also provide a level2 implementation, pick a bailing condition that works fine on my laptop but no further investigation on other machines. We do not list all the level2 implementations below, but they are a part of this PR.
{c,z}hetrrk: A complex hermitian triangular rank k update assuming the update matrix is triangular (thank you Mathieu Faverge for the naming!){s,d}sytrrk: A real symmetric triangular rank k update assuming the update matrix is triangular{s,d,c,z}trmmoop: An out of place triangular-general matrix product. Similar to a trmm + an add.{s,d,c,z]trmvoop: An out of place triangular-vector product. Similar to a trmv + axpy{c,z}acxpy: computing{c,z,s,d}trsm_mod: A modified triangular solve that allows for solving singular systems with the assumption that A has a zero row/column where there are explicit 0s on the diagonal and that B has the necessary values equal to 0. We do not enforce this condition in the code. The purpose of this is to be able to use the newlarft_utscheme by applyinglarft_uton entry withapply_t = 's'{c,z,s,d}trtri_mod: A modified recursive triangular inverse to mimic the behavior oflarftwhen inverting the triangular factor before returning.{c,z,s,d}trti2_mod: A modified level2 triangular inverse to mimic the behavior oflarft_lvl2 when inverting the triangular factor before returning.The previous two routines are necessary when we generate a reflector that doesn't do anything (IE:$\tau=0$ from
larfg) There are two main cases we found that use this:or/ungqr)When the triangular matrix is non-singular (ie all$\tau$ values are non-zero), then it is equivalent to the existing routines
trtriandtrti2Another option would be to change$\tau=2$ , but this would cause a triangular input to have it's diagonals sign's flipped, which is not necessarily desirable.
larfgto return a reflector withThe following existing routines were modified as well:
{c,z,s,d}larft:LQandRQfactorizations. These paths are not yet needed, but will be a part of a future PR where we optimize the Q generation family, I felt it best to include these paths here as it affects the triangular factor generation. I can wait to add these paths until we explicitly need them if desired.larftto uselarft_utas some experiments we have show the ut-based version is more performant for certain matrix sizes. This can be removed or changed back to level2 if desired.Checklist