Repository navigation
xORMQR and xLARFB discrepancy with documentation #1011
Description
Activity
For LARFB, you are correct. We could change: "The elements equal to 1 are not stored; the corresponding array elements are modified but restored on exit. The rest of the array is not used." to "The upper triangular part of V (including its diagonal) is not referenced." Thanks. Feel welcome to provide a PR for LARFB.
For ORMQR, you are correct too. If we use the "unblocked" version of the code, we use ORM2R and then "the diagonal elements of A are modified but restored on exit". And then after line 96 in ORMQR write "A is modified by the routine but restored on exit."
Also, in ORM2R and ORMQR, A could be labelled as "input/output" and not "input". Well, I am not sure about this. A is the same in input and in output, but it is changed and that would help for thread safety.
However we should really "fix" ORM2R (by fixing LARF) so that LARF does not require the "1" on the diagonal of V so that ORM2R does not "modify and restore on exit" entries of V. Maybe we create a new routine for this, such as LARF1. Or just change LARF proper.
Another fix for ORMQR would be to remove the call to ORM2R and always call the "blocked code". So remove
IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN * * Use unblocked code * CALL DORM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, $ IINFO ) ELSE * * Use blocked code * ***** [ KEEP THIS PART ] ***** END IF
So that we never use ORM2R in ORMQR. That would be do the trick.
If you want to propose a PR for the documentation of ORMQR, this is OK. My preference is to fix LARF, I can work on this but this will have to wait for a month or so.
Reacted by Eduard Fedorenkov- added a commit that references this issue
on May 21, 2024 @langou, about ORM2R and LARF. Looks like the simple solution is this:
- ORM2R: double workspace for ORM2R. One part is for copy of V and another for storing C**T * V. Finaly remove code which modifies A.
- LARF: just add COPY in LARF to copy V inside work with something like
work(1) = 1;.
Also, I think that removing of unblocked algorithm ORM2R from ORMQR is not so bad idea. But, may be there are some reasons that can damaged some tests. I don't know.
about ORM2R and LARF. Looks like the simple solution is this:
(1) ORM2R: double workspace for ORM2R. One part is for copy of V and another for storing C**T * V. Finaly remove code which modifies A.
(2) LARF: just add COPY in LARF to copy V inside work with something like work(1) = 1;.Thanks for sharing ideas. Some comments: () Increasing necessary workspace is not a good idea in general and we have been reluctant to do it for backward compatibility. () ORM2R will be fixed if we fix LARF. (*) For LARF, we should write a version that assumes that the first or the last component of v is a 1. LARF1F (1 first) and LARF1L (1 last) for example. And then propagate to ORM2R. I will try to get this done by this Summer.
Also, I think that removing of unblocked algorithm ORM2R from ORMQR is not so bad idea. But, may be there are some reasons that can damaged some tests. I don't know.
This is a good idea but I am leery to go in this direction. This kind of change is up to people tuning the code. I think the code as is fine. The condition to call ORM2R (inside ORMQR) is
Line 275 in 4a8ed6b
and I agree that this can be changed. This should not damage any test. (As long as this is done correctly.) I agree that changing this line of code could bring some speedup. That being said, I do not think we should change this for now.Julien.
Reacted by Eduard Fedorenkov@jprhyne just pushed his implementation in #1020
I am sorry for the duplication of work. I did not expect it.
It seems like both implementations are suffering from corner cases.
While we need to get these corner cases working, I do not want to create stress and start a race on who finishes first. It would be great if both of you are able contribute to LAPACK through your work. Maybe we can agree that you will work on LARF1F, and the other work on LARF1L. Or someone does C and S and the other one does Z and D. If that's acceptable. I think the split {C,S} for one and {Z,D} for the other is the best for now.
The two codes (@jprhyne's #1020 and @EduardFedorenkov #1019) look similar-ish that I feel it is possible to merge both work and it would be great to have both of you contributing to the software.
In any case, it is great to have contributions and a start of resolve for the issue. Thanks to both of you for the work.
Reacted by Eduard FedorenkovHi, @langou and @jprhyne!
It's OK, let's work together!Great idea. I support your suggestions! I will try to find the solution for corner cases as well as @jprhyne
let's split {Z,D} and {C,S}. Also, we need to align code. So, @jprhyne can review my PR and I can review #1020. When we all agree about code we can just merge #1019 with {Z,D} and #1020 with {C,S} for example.Related, not necessarily urgent, but it would be good to push the change to more subroutines than only ORM2R. For example GEQR2 as the same "guilty" lines:
Lines 184 to 191 in dd2e5ef
Let us fix ORM2R first.
I did a grep of LARF and found all these routines:
xGEBD2
xGEHD2
xGELQ2
xGEQL2
xGEQP3RK
xGEQR2
xGEQR2P
xGERQ2
xLAQP2
xLAQP2RK
xLAQR2
xLAQR3
xLARFX
xOPMTR
xORBDB
xORBDB1
xORBDB2
xORBDB3
xORBDB4
xORG2L
xORG2R
xORGL2
xORGR2
xORM2L
xORM2R
xORML2
xORMR2Note: the xOR prefixes are for real arithmetic {S,D} and so the corresponding complex arithmetic subroutines would be under xUN for complex arithmetic {C,Z}.
Reacted by Eduard Fedorenkov- added 4 commits that reference this issue
on May 31, 2024 12 remaining items
I think this PR for {S, C} is about 95% done. I look over one more time and add some small fixes. I have some notes to discus:
- We talk with @jprhyne about zero vector as input for xLARF1F and decide that the smallest LASTV must be limited by 1. But unfortunately I have no idea how to deal with this issue in xLARF1L (I mean zero vector input)
- Interesting case is in xLARFX. Looks like changes in xLAEXC is needed in this case. And maybe xLARFX must be spited by xLARFX1F and xLARFX1L as well
- added 4 commits that reference this issue
on Jun 7, 2024 Related to @EduardFedorenkov point about lastv with larf1l. I have 2 questions/ideas about how to proceed.
We could
1: Since we are assuming v(lastv)=1 (in fact we guaranteed this was the case in the previous paradigm) we could replace it with a firstv parameter and start scanning at the start of v.
2: scrap the larf1l.f files and instead call larf1f.f with a negative incv (I have not really tested this so I am not too sure how feasible this is). This has the benefit of one less extra file but could lead to some headaches depending on if there are edge cases that don't play nice with this paradigm. My naive attempts of just dropping in dlarf1f where dlarf1l was and negating incv so far have added some failures but it may just need us to think differently about the incv < 0 case.edit: fixed incorrect file names
After some prototyping, I don't think option 2 is the best idea. I'm not too familiar with how the negative incv values should behave, and the ideas that I can think of to use this are harder to read than the current implementation so I'll be not working on that anymore, but maybe someone can look at it at some point and revisit this idea.
Until then I'll be working on the 1l with a "firstv" variable instead of a lastv
edit: Fixed formatting accidentally referencing another post
the 1l with a "firstv" variable
I had not thought about this. Oh my.
@jprhyne is making a good point and proposing to introduce a new feature that is not in LARF and that is an advantage of this LARF1L / LARF1F. The current LARF code always scan and look for zeros from the bottom (last entry of the vector). However when the "one" is at the bottom, it makes more sense to scan from start of the vector. What @jprhyne is proposing would be a new feature. (If I understand correctly.) And this makes sense.
I think that the better way is to change the set of the arguments in LARF1F and LARF1L.
What about to add parameter that store index in which ONE is assumed to be. No need to compute lastv or firstv in this case. What do you think @jprhyne, @langou?
Because even the current realization of LARF1F failed the tests if in place where we expect ONE will be the pure ZERO.
I like this idea, but I don't see how it will be different for larf1f
The current idea behind my attempt is to assume v(1) = 1 without referencing it. Something isn't playing nice with that scheme and is probably a simple indexing issue with the lastv scanning or something else that I need to just sit down and go through slowly with many cases, so adding another parameter shouldn't change anything about the algorithm itself. We would just do
DO WHILE(LASTV.GT.FIRSTV_ARGUMENT .AND. V(I) = ZERO) LASTV = LASTV-1 I = I - INCV END DOor maybe
DO WHILE(LASTV.GE.FIRSTV_ARGUMENT+1 .AND. V(I) = ZERO) LASTV = LASTV-1 I = I - INCV END DOAnd later on we would instead of referencing
V(1)andC(1,1)we would have
V(1 + (FIRSTV_ARGUMENT-1) * INCV)andC(FIRSTV_ARGUMENT,1)(and the C parameters flipped when SIDE is changed)Which is exactly what the form of the
FIRSTVscanning I was discussing above would do. (Which is when I noticed a bug I pushed a small fix for earlier by removing this scanning entirely before continuing to work on it). I'm really not sure how helpful this would be as personally these scanning parameters are not intuitive and our tests are designed around this behavior not existing. So I will be trying to get everything else done before that functionality.Then this raises another question of since we're adding another parameter why not have a character input
STOREV(like how is done inxLARFT) to merge these together. I'm not sure if that's the best of ideas as at this point we're makingxLARF1{F,L}somewhat more complicated than it needs to be for applying a single householder reflector with a small change in structure. But, I'm not sure what the correct move is for this.- added 2 commits that reference this issue
on Jun 11, 2024
Let's suppose that we have a task to apply vectors from A (returned by xGEQRF) to the set of matrices by xORMQR in parallel. I found that current LAPACK has discrepancy with documentation that introduce doubts when you try to solve this task.
Maybe it's worth doing some corrections in LAPACK documentations for xLARFB and xORMQR?