[Request for comments] Properly solving the *TRTRS singularity detection problem about lapack HOT 2 OPEN

option # 1 is for all intentions and purposes, impossible. That ship has sailed a few decades ago.

Some users could end up faced with runtime failures ( INFO > 0) after such hypothetical changes to LAPACK, although it could be argued that it is better that they become aware that their matrices are potentially numerically unstable.

This one is easy. That's a clean no (we did this in SciPy it wasn't fun). They survived with potential singular arrays all the way here. You are just going to annoy the whole world about a no-op warning and get a bad rep for no particularly strong reasoning.

Secondly, it would still not let users customize/adjust the threshold where an error is signaled. Contrast this with *GELSY: it does pivoted QR and asks the user to specify RCOND, where RCOND is used to determine the effective rank of A.

Correct, and one of the biggest complains about this, is where you get the rcond from. So it won't solve your problems, in fact you now have to guess a rcond for your array or perform another painful long long command since ?##con wants the ##trf result to compute and the function signatures are unnerving just to solve a single AX = B equation.

For # 2

I always wanted to have this anyways, so don't get me wrong. But the question is how are you going to do it? Just by checking the absolute value is not enough. So you should get the condition number for inversion. Typical application of this is you do one ?##trf, then squeeze in a ?##con, and then continue with ?##trs if con result goes OK (Or I don't know any other way to implement it). But note that rcond computation is not exact and not really that fast(calls latrs); so you are approximating it. Hence there will be some art to it and hence false positives.

In my opinion, if there is any chance to touch these algorithms, it has to implement the matlab backslash behavior. Check out the flowchart here. This is not really a complicated algorithm and allows for quick for loops over the array.

That is to say, it discovers the structure of the array by running down the entries, checking for bi, tri diagonal, upper/lower triangular, symmetric and hermitian, and pos def. and decide the algorithm and the ?xxcon check as it discovers choosing the easy way out.

The rest is, it seems to me, just unnecessary labor for the users, legacy of 40 year old habits with 5 different calls with respective lworks and aligning memory for the results and so on. It should get easier not more laborious by the decade.

That would be my door # 3

from lapack.

option # 1 is for all intentions and purposes, impossible. That ship has sailed a few decades ago.
You are just going to annoy the whole world about a no-op warning and get a bad rep for no particularly strong reasoning.

I thought so as well, and initially rejected the idea myself as well. But still, I am curious to know the community attitude is towards the idea. I mean, who knows, some communities are lot more gung-ho about minor changes to function behaviour, so might as well put it on the table.
But I agree, my first instinct is also to avoid option # 1.

Secondly, it would still not let users customize/adjust the threshold where an error is signaled. Contrast this with *GELSY: it does pivoted QR and asks the user to specify RCOND, where RCOND is used to determine the effective rank of A.

Correct, and one of the biggest complains about this, is where you get the rcond from. So it won't solve your problems, in fact you now have to guess a rcond for your array or perform another painful long long command since ?##con wants the ##trf result to compute and the function signatures are unnerving just to solve a single AX = B equation.

I used *GELSY just as a precedent for having some kind of decision theshold, I 100% agree that having users sit there scratching their head, trying to figure out where to get RCOND from to call *GELSY, is a frustrating experience.
But I think the biggest reason for that is the lack of a sane default. This is why I suggested putting in a default that ought to have very few false positives.

But the question is how are you going to do it? Just by checking the absolute value is not enough. So you should get the condition number for inversion.

I am not convinced that going through the full condition number estimation is necessary for triangular solvers like *TRTRS. For a triangular matrix the eigenvalues are just the diagonal elements, and the determinant is just the product of the diagonal elements. And there is a convenient "natural" default to check against in the form of subnormality. Beyond the individual elements one could also check the determinant as it is so easy to compute for triangular matrices. If any of the diagonal elements or the determinant is subnormal, trigger the error exit.

from lapack.