[C/Z]LAUUM / [C/Z]LAUU2 only reads the real part of the complex values on the diagonal of input matrix

[langou]

This is related to work of @jprhyne with #1402
And discussion from @ACSimon33 at #1373 (comment) to use LAUUM in LARFT

So we are reviewing [C/Z]LAUUM / [C/Z]LAUU2, it appeared that [C/Z]LAUUM / [C/Z]LAUU2 only reads the real part of the complex values on the diagonal of input matrix

See line 164 at:

lapack/SRC/zlauu2.f

Lines 163 to 178 in 8835bfa

	DO 10 I = 1, N
	AII = DBLE( A( I, I ) )
	IF( I.LT.N ) THEN
	A( I, I ) = AII*AII + DBLE( ZDOTC( N-I, A( I, I+1 ),
	$ LDA,
	$ A( I, I+1 ), LDA ) )
	CALL ZLACGV( N-I, A( I, I+1 ), LDA )
	CALL ZGEMV( 'No transpose', I-1, N-I, ONE, A( 1,
	$ I+1 ),
	$ LDA, A( I, I+1 ), LDA, DCMPLX( AII ),
	$ A( 1, I ), 1 )
	CALL ZLACGV( N-I, A( I, I+1 ), LDA )
	ELSE
	CALL ZDSCAL( I, AII, A( 1, I ), 1 )
	END IF
	10 CONTINUE

And so, @jprhyne and I were wondering whether that was a good idea or not. This "feature" is not even described in the comment of the routine. (@jprhyne plans on fixing this by adding some comments if we leave it like that.)

On the hand, one application of LAUUM is multiplication of Cholesky factors and so the diagonal is real. So it does not hurt to set the imaginary part of the complex number to zero. In our LARFT application, the diagonal of the triangular matrix is also real since the triangular matrix is unit. So there is no real harm to leave it as is.

But this is weird. It feels like some users might want to form U U^T where U has a non-zero imaginary part for the diagonal elements. I think it would not be too hard to do and give the routine a broader useage case.

I think if there is no broad support, @jprhyne and I plan to leave the feature in and document it.

If support for changing things, let it known.

[jprhyne]

One small difference that I noticed is the routine also requires the diagonal explicitly have 0s in the imaginary component to behave in a way compatible with the Cholesky factors. I'm not familiar with the potrf family, do they explicitly set the imaginary component to 0 before returning or is it implicit?

If we set $n=1$ then this is just a wrapper around ZDSCAL and computes A = REAL(A) * A so whatever is in the imaginary component also gets scaled. I think that we should do either

1. not require the imaginary component be explicitly set to 0 (or do it ourselves since it will be overwritten anyway)
2. don't make any assumptions on the data of the diagonal of A

[ilayn]

I think if a user wants a triangular multiplication they should use trmm and should not repurpose this. Hence probably these need a documentation update and not a behavior change.

[jprhyne]

I can see that if we update the docs, we might as well just zero out the imaginary component of the diagonal elements. This ensures that when we call other routines like zdscal (or ztrmm inside of zlauum.f) we don't get inconsistent behavior if the user modified the imaginary component of the diagonals.

I'm only in favor of adding zeroing the imaginary components of the diagonals since (1) we are already overwriting this value and (2) we delegate to routines that require this zero to be explicitly set.

It is an extra pass over the diagonal of the data, but I omitting that level of optimization is worth the extra explicitness on the desired behavior.

Since we are only grabbing the real component of the diagonal I would expect this behavior (ignoring whatever is in the imaginary) to be mimicked everywhere viable or set explicitly if needed (akin to calling LACGV before the GEMV in LAUU2)

[langou]

@jprhyne and I reviewed the current code in LAPACK for LAUUM and LAUU2 and this is not good. We need to change it.

In some places, the current code zeroes out the imaginary part of diagonal elements. In other places, the code works on the diagonal elements "as is".

As a consequence, the only way to have a meaningful results is to call with diagonal elements with zero imaginary parts.

Having with diagonal elements with zero imaginary parts is the only use case that we have encountered so far and this is why the routines have been operating as such.

If you call with nonzero imaginary parts on the diagonal, the routine is currently returning nonsensical results. It does not zero out for you. (It does at some places but not all.) Since it sometimes zeros out, sometimes not, the result is nonsensical.

We could say that the only correct way to call LAUUM/LAUU2 is with zero imaginary part on the diagonal but all this does not make much sense. Since it is easy to enable a sensible behavior.

What would make sense is either consistently zeroing out the the imaginary part of diagonal elements everywhere in the code (that would be akin to POTRF or HERK). Or zeroing out imaginary part on the diagonal nowhere on the code (that would be akin to TRMM, TRSM, etc.).

Either would not disrupt standard operation for matrices who have real diagonals.

I do not think the zeroing out of imaginary part of diagonal elements make any sense for this application. This is done in BLAS/LAPACK when we assume that we are dealing with the upper part (or lower) of a Hermitian matrix (e.g. HERK, POTRF). We can see having this behavior for the Cholesky factor, but (1) a Cholesky factor should have the imaginary part of diagonal elements zeroed out, and (2) it is not a good idea to specialize a routine so much (in my opinion), if it can as well handle general triangular matrices.

My preference would be to remove the automatic zeroing out of the imaginary part of diagonal elements in these few places in the code so that we have a nice "triangular HERK in place" algorithm.

Thanks for point this out @jprhyne. Can you please fix it?

[jprhyne]

Yes I can, I'll incorporate it into PR #1402.

[jprhyne]

To be clear, I support not zeroing out the diagonal elements before operations, but ensuring they are exactly real on exit is I think good since that is compatible with HERK