xORMQR and xLARFB discrepancy with documentation

[EduardFedorenkov]

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.

1. xORMQR : According to documentation the unblocked code branch - xORM2R modify diagonal of A matrix and restore it on the exit. But xORMQR has no such remark about A matrix. So, in this case we can't use xORMQR in parallel because of rewriting the diagonal of A.
2. xLARFB : According to documentation vectors in V are modified and restore on the exit. But according to the code GEMM and TRMM with diag='Unit' don't touch the diagonal elements of V.
   

Maybe it's worth doing some corrections in LAPACK documentations for xLARFB and xORMQR?

[langou]

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.

[EduardFedorenkov]

@langou, 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;.

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.

[langou]

Hi @EduardFedorenkov,

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

lapack/SRC/dormqr.f

Line 275 in 4a8ed6b

IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN

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.

[langou]

@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.

[EduardFedorenkov]

Hi, @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.

[langou]

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:

lapack/SRC/dgeqr2.f

Lines 184 to 191 in dd2e5ef

	*
	* Apply H(i) to A(i:m,i+1:n) from the left
	*
	AII = A( I, I )
	A( I, I ) = ONE
	CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ),
	$ A( I, I+1 ), LDA, WORK )
	A( I, I ) = AII

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
xORMR2

Note: 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}.