NaN check in LAPACKE_?lantr isn't correct for m != n

[ACSimon33]

Description

The NaN check in LAPACKE_?lantr isn't checking the entire matrix for non-square matrices. The check is implemented like this:

if( LAPACKE_ztr_nancheck( matrix_layout, uplo, diag, MIN(m,n), a, lda ) ) {
    return -7;
}

There are two solutions to this. Either we add another check for the part which isn't checked - that could look like this for example:

if( LAPACKE_ztr_nancheck( matrix_layout, uplo, diag, MIN(m,n), a, lda ) ) {
    return -7;
}
if( LAPACKE_lsame( uplo, 'u' ) && n > m) {
    lapack_int offset = m * ((matrix_layout == LAPACK_COL_MAJOR) ? lda : 1);
    if( LAPACKE_zge_nancheck( matrix_layout, m, n - m, &a[offset], lda ) ) {
        return -7;
    }
}
if( LAPACKE_lsame( uplo, 'l' ) && m > n) {
    lapack_int offset = n * ((matrix_layout == LAPACK_ROW_MAJOR) ? lda : 1);
    if( LAPACKE_zge_nancheck( matrix_layout, m - n, n, &a[offset], lda ) ) {
        return -7;
    }
}

The cleaner solution in my opinion would be to change LAPACKE_?tr_nancheck to accept non-square matrices. That should be fairly simple by just adding a second size parameter and adjusting the loops. The NaN check in LAPACKE_zlarfb would benefit from this as well since there we have the same problem (which is already addressed via the first strategy).

The only downside of changing LAPACKE_?tr_nancheck is that there would be a lot of files that need fixes due to the interface change (64 files).

I'd like to know which strategy you would prefer.

Checklist

• I'd be willing to make a PR to solve this issue

[langou]

Hi,

This is a good topic of conversation. I am not sure what a nonsquare triangular matrix is. If we think about a rectangular m-by-n matrix and we say it is triangular, I feel it is not obvious where the triangle is. Is it above or below the rectangle, left or right of the rectangle? I feel different people will give different answers. I think if we look at LQ, QL, RQ and QR, we will have all these shapes naturally appearing.

I do not think we have nonsquare triangular matrices in LAPACK or in BLAS. I think triangular matrices are always square in LAPACK and BLAS.

Although I think that I see what you mean with xLARFB. Yes the householder reflector vector matrix is m-by-n (m>=n) lower triangular unit. So to check NaN in it we need to check NaN in the strictly lower triangular and then the rectangle below (if it exists, case m>n+1).

For now my preference would be to keep LAPACKE_ztr_nancheck as it is, and then to append a zge check. I can be convinced otherwise though. Happy to change my current opinion.

Julien

[ACSimon33]

You're right, that the term triangular matrix always refers to a square matrix. I wasn't aware of that. The term I meant to use is trapezoidal matrix. This kind of matrix is also mentioned in the LAPACK documentation of some routines, e.g. xLANTR, xTZRZF, xLATRZ, xGETRF, etc. They all use the same definition of trapezoidal where the triangle is left or above the rectangle (see the link).
There are places where the definition is different, e.g. for xLARFBwith DIRECT = 'B', but there the documentation doesn't refer to the matrix as trapezoidal.

We could also introduce a new NaN check for trapezoidal matrices (e.g. LAPACKE_ztz_nancheck) and keep LAPACKE_ztr_nancheck as it is. But I think I'll create a PR with the first strategy for now, since I already implemented that.

[langou]

Hi @ACSimon33,

Thanks for bringing up xLANTR. This is a good point. I wish xLANTR was named xLANTZ.

There is a difference between a norm computation and a NaN check. For the 1-norm computation, you cannot reduce the two output of xLANGE and (proper) xLANTR to get the output of xLANTZ. (You can easily do so for the Frob, inf, 1 and M norms.) Whereas for nancheck, you always can easily reduce. So I see why xLANTR is handling some trapezoidal matrices, but I wish it would be called xLANTZ.

I liked the Z in the naming of xTZRZF, xLATRZ.

A trapezoidal NaN check would be useful. Thanks for bringing this up.

My preference is (second strategy) introduce a new nancheck subroutine (e.g. LAPACKE_ztz_nancheck) and keep LAPACKE_ztr_nancheck.

That being said, (1) I am thankful that you are considering doing a PR with trapezoidal NaN check, (2), as you pointed, there is a precedent with xLANTR, (3) you are the one who spent time on it, and are designing the PR, so I am very comfortable if we end up with (first strategy) LAPACKE_ztr_nancheck that checks trapezoidal matrices. That's 100% A-OK with me.

Thanks!