ralna/spral

Bug in ssids_enquire

Closed this issue · 0 comments

I have been using the ssids_enquire_indef function to provide the (block) diagonal entries, d, of the factors. I frequently find that the reported values are infinite. On looking closely at the relevant c++ code NumericSubtree.hxx, line 454, I find that this is where the data is extracted from the internal structure. There is an if-else that starts with if(i+1==nelim || std::isfinite(dptr[2i+2])). In the case where this test fails (i.e., the else), then dptr[2i+2] is infinite ... but the third (d++) statement sets d++ to this infinity. I frequently see these infinite components in my output. My only supposition is that this should be (d++) = dptr[2i+3] not 2i+2; without digging further, I can only guess that, in each column, the diagonal and subdiagonal is stored consecutively, that the 2nd column in a 2x2 block is flagged by having an inf as its first value, and that the other (symmetric) value occurs as the second value. But then that is only a guess, and there are no comments to help! If I set line 454 to
(d++) = dptr[2i+3]; /* not 2*i+2 as stated ?? */
I get the finite "right answer", but whether this is a fluke is anyone's guess. I have run the "make check" tests and all pass, but they did before, and I suspect the tests are not very comprehensive with these peripheral utilities.