BLAS and LAPACK Subroutines of Real Skew-symmetric Matrix

[sh-zheng]

This PR includes some new subroutines of real skew-symmetric matrix. These subroutines include BLAS2/3 operations, linear solver and eigensolver in LAPACK.

Please refer to the attachment document for more details.

lawn_skew_symmetric.pdf

[sh-zheng]

Hi~, are there any comments, or further things to do about this PR?

[sh-zheng]

Hi, @mgates3 , I have made some updates, including the dates, authorship and the doxygen doc. Could you please continue the review?

[mgates3]

To be honest, this PR is too large to review. It needs to be broken into smaller, manageable chunks. Even then, with such a large change, it should be coordinated with LAPACK maintainers beforehand, e.g., to discuss naming conventions. I'll discuss this PR with other LAPACK maintainers at our next meeting.

[langou]

I agree with @mgates3.

I agree that it would be better and easier too and probably faster to break down this PR in several smaller PRs that we would slowly integrate. I think the overall work is great and it is important to give context to the contribution and this PR gives plenty of context for each of the contribution. But going with a finer resolution with the integration could be better. (And not necessarily slower.)

For example, maybe we can start by KYMM.

If we go that road, we can continue refer to this PR as needed. So this PR is not lost and is important to enable a way forward.

With respect to KYMM, I have been thinking whether we want it in the BLAS or in LAPACK. After the Level 3 BLAS paper (1990), BLAS-like routines have been integrated in LAPACK. Only recently (#887) were new interfaces added to BLAS. All in all, I am leaning to have KYMM part of BLAS, but it would be good to hear from others.

With respect to KYMM, I wonder if KEMM is not better than KYMM. I assume the Y is for sYmmetric. But an E for hErmitian might be more appropriate.

One though about this PR is that this is real arithmetic only and there is no support for complex arithmetic. I do not know what to do about this thought.

In the meantime, @sh-zheng, if you are willing to take a more incremental approach, maybe a parallel PR with only KYMM So BLAS/SRC, BLAS/TESTING, CBLAS, DOCS, then maybe this would be one way to make progress. This would probably be less intimidating for all.

An incremental approach would be more manageable. This also might result in a faster turn around time for the whole PR in the end as eluded to by @mgates3.

Comments welcome.

[sh-zheng]

Thank you so much for the review @mgates3 @langou .

I agree with that this pr should be divided into finer sub-prs to improve the clarity. I plan to arrange the sub-pr in following 5 stages.
a) Skew-symmetric blas and its testcases.
b) Corresponding cblas.
c) Skew-symmetric lapack linear-solver and its testcases.
d) Skew-symmetric lapack eigen-solver and its testcases.
e) Corresponding lapacke api.
And this pr can be kept as a code base of the whole large submission.

About the belonging of kymm (also kymv, etc.), I tend to put them in blas, for the functional duality of api (symm vs. kymm). The different distribution (symm in blas and kymm in lapack) may lead to potential confusion.

About the naming rule, I tend to seperate the real skew-symmetric ("ky") and the skew-Hermitian subroutines ("ke"), which are dual with real symmetric ("sy") and Hermitian subroutines ("he"). I haven't started to realize the skew-Hermitian code because the real skew-symmetric code may be more important in application, and maybe most of skew-Hermitian demand can be implemented by multiplying i.

[ilayn]

Do we still need to encode into two letters a lot of meaning? Why not just (s,d,c,z)skewsymm ? Or even plain (s,d,c,z)skew_hermitian for that matter?

[sh-zheng]

Do we still need to encode into two letters a lot of meaning? Why not just (s,d,c,z)skewsymm ? Or even plain (s,d,c,z)skew_hermitian for that matter?

Although the trend of naming rule is more loosed today, for such a "dual" subroutine (eg. symmetric vs. skew-symmetric), a traditional name with obvious contrast may be more clear for users.

[mgates3]

skewsymm and skewhemm make a lot of sense to me — very clear (assuming you already know symm and hemm), whereas kymm and kemm are pretty cryptic (sKew-sYmmetric, sKew-hErmitian).

[sh-zheng]

skewsymm and skewhemm make a lot of sense to me — very clear (assuming you already know symm and hemm), whereas kymm and kemm are pretty cryptic (sKew-sYmmetric, sKew-hErmitian).

Since the naming rules are important to users, it may need wider discussion at meetings. If some majority opinion has been made, I would change the naming rules in the following separated pr.

[ilayn]

That is true but kymm/kemm are also not decided by users. Hence there is always some decision priority by the code authors regardless. Moreover users typically don't have any preference until something is inconvenient or unergonomic.

[sh-zheng]

I reviewed the naming scheme of the whole skew-symmetric submission, and still tend to believe that the traditional naming style "ky" is a suitable choice (at least in the present). This tendency comes from the following considerations.

a. Users' habits mean that they are unlikely to be confused by this naming scheme. A survey in lawn 290 showed that most users referred to lapack user guide (including document of mkl, cublas, etc.), and the mapping between matrix types and prefix of name is one of the most prominent parts in any document. So a user who frequently uses certain functions might soon get used to this naming scheme. That is to say, users may not even pay attention to this naming scheme again in their subsequent use.

b. Although the naming scheme of blas/lapack came from the fixed format fortran, and in the free format fortran it was relaxed, the prefix (data type and matrix type) of the name of subroutines submitted later still obeyed the traditional scheme, and just relaxed the part of specific functions. For example, ssysv, ssysv_rook, ssysv_rk, ssysv_aa, etc. In this submission, if the name of skew-symmetric subroutines are expanded into completely plain form (like skewsymm), it may result in two potential problems. Firstly, users' habits and preconceived notions will be broken and they may need to adapt to a totally new format. Secondly, some internal code that parses calls based on subroutine names in blas/lapack will no longer be applicable, and it has to be adjusted and tested to adapt the completely plain naming scheme.

c. Finally I would like to explain the choice of "ky". In pfapack library the prefix "sk" was adopted. But our choice may be more suitable in reserving namespace for tridiagonal, band, and packed storage subroutines. So in this submission the "k" was extracted as the most fundamental feature of matrix type. Then another possibility "kw" emerges ("ke" is not considered due to the confusion with the Hermitian case). However if we take the duality between the symmetric and skew-symmetric cases into consideration, "y" should not be discarded. This is because in shared code, especially environment configuration and testcases, the symmetric and skew-symmetric code can be reused and we just need to swapping the distinguishing letters between "k" and "s". What's more, the correspondence between st vs kt, sb vs kb and sp vs kp, may be more natural in semantic level.

[ilayn]

I don't think users had any choice in the matter for a very long time as these names were forced upon them by technology. Plus there were very little choices to be made ge, sy, po, gb were quite distinctive and intuitive since po was positive sy was symmetric.

However note that this was still a hack to bypass a naming hurdle. Once the hurdle is removed all there is left was habits. This is a very difficult thing to overcome but there is no reason to continue technically. I think half of the LAPACK algorithms are already badly named just because of this habit. Everytime I have to check whether it is lagrf or lagrp, or lagrfp, or whatever in the same day just to find the right one. With no offense ky, ke is not helping either.

We have no reason to continue this naming and I need to be convinced that there are users who wants to decrypt a puzzle instead of reading the name of the routine directly. I

In fact, I am positive that this naming scheme is a self-inflicted pain as soon as the file name limitations are removed in the 90s.

You might not like the name suggestions and that's perfectly fine but I think we can just lay to rest the argument that users will want the old scheme. I don't think we have any data supporting or against it.

Firstly, users' habits and preconceived notions will be broken and they may need to adapt to a totally new format.

Based on the argument above I don't think so. Aasen or rook did not break anyone's expectation.

Secondly, some internal code that parses calls based on subroutine names in blas/lapack will no longer be applicable, and it has to be adjusted and tested to adapt the completely plain naming scheme.

That's a very old Fortran77ism that we shoupd not be using anyways. If you need an NB parameter you can always base it on existinf routines or even better use an explicit number independent from NB.

[mgates3]

Regarding skew-symmetric vs. skew-Hermitian, ky and ke could work for skew-symmetric and skew-Hermitian, but kb, kp, kt would conflict between skew-symmetric and skew-Hermitian:

sp and hp => kp
sb and hb => kb
st and ht => kt

Whereas adding skew would avoid that conflict (skewsy, skewhe, skewsp, skewhp, skewsb, skewhb, etc.).

Instead of ky, we could use ke to mean skew-Hermitian, and consider real-skew-symmetric as the same as real-skew-Hermitian. Then ke, kp, kb, kt would all mean skew-Hermitian, with the understanding that in real, skew-Hermitian == skew-symmetric. For instance, for Cholesky in SLATE we use herk for the update, with the understanding that in real, herk means the same as syrk. (This overloading is essential for C++ templating.) But this would not provide a name for complex skew-symmetric matrices, which an Internet search does yield some hits for.

Yes, ilaenv would need to be updated for skew matrices, no matter what naming scheme is used, but this would be easy. Are there other codes or scripts that assume letters 2-3 of the name are the matrix type?

There are some routines already that encode the matrix type differently, not in characters 2-3, such as dsposv, dsgesv, lan{ge,gb,he,hf,ht,tr,...}, and la{he,sy}f*. (Though I also think some of these aren't very good names. {ge,...}norm would be much nicer.)