bugGNU Octave - Bugs: bug #63785, float compared to 0 in __ichol__.cc


bug #63785: float compared to 0 in __ichol__.cc

Submitter:  Dmitri A. Sergatskov <dasergatskov>
Submitted:  Fri 10 Feb 2023 09:10:46 PM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 8.0.90 Operating System:  * Any
Fixed Release:  8.3.0 Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Wed 19 Jul 2023 04:53:19 PM UTC, comment #7: 

I checked in the changeset here: http://hg.savannah.gnu.org/hgweb/octave/rev/b73ebf017461.

My feeling is that we are working around CPU vendors who should be doing a better job of maintaining consistency of results regardless of how a high-level instruction like multiply is mapped onto individual hardware units.

If the errors do grow beyond eps, the threshold in this changeset, then we will start to get new reports and will have to develop a different solution.

Marking as Fixed and closing report.

Rik <rik5>
Group administrator
Tue 18 Jul 2023 09:48:20 PM UTC, comment #6: 

I guess we can close it now.

I am still uneasy about it. With sse the test never fails, but with AVX one can construct the matrix (scale it by large number and increase the size) that will fail. It feels that we hide the bug rather than fix it.

I was trying to construct a standalone C++ example that demonstrate the problem to ask compiler people about it, but I was not able to.


Dmitri A. Sergatskov <dasergatskov>
Tue 18 Jul 2023 09:37:13 PM UTC, comment #5: 

I took a look at the magnitude of the imaginary term and it is very small (1e-19).  Hence, I propose to compare to eps() rather than scaling the tolerance by the size of the real pivot.  Patch attached for review.

(file #54938)

Rik <rik5>
Group administrator
Wed 08 Mar 2023 08:32:31 PM UTC, comment #4: 

Oops I forgot that k6-8350 already has avx (just no avx2).
have to make sure no AVX optimization is enabled.
So nothing new.


Dmitri A. Sergatskov <dasergatskov>
Wed 08 Mar 2023 07:26:37 PM UTC, comment #3: 

With gcc 13.0.1 in fedora
this (and pcg) tests are failing with default optimization gcc now
on amd k6-8350 (fedora buildbots).


Dmitri A. Sergatskov <dasergatskov>
Sun 12 Feb 2023 08:23:53 PM UTC, comment #2: 

May be this is a gcc bug. I played with code:

for ii = 1:100 a = randn(N) + i*randn(N); b = R*sparse(a'*a); ii, ichol(b); endfor

(with different N and R)

and situation gets progressively worth for the larger N.

In clang compiled code (with the same optimizations)
Im part is strictly 0.

E.g. gcc (N =4, R=1):

ii = 29
Real: 5.37065 Im: 0
Real: 5.05964 Im: 8.51589e-19
Real: 7.94997 Im: 5.10252e-17
Real: 0.0602546 Im: 4.3797e-16
error: ichol: non-real pivot encountered.  The matrix must be Hermitian.
error: called from
    ichol at line 246 column 9


Dmitri A. Sergatskov <dasergatskov>
Sun 12 Feb 2023 07:06:10 PM UTC, comment #1: 

I think in general you are right that we shouldn't be comparing exact values (derived from pure math) with calculated values made on computers with finite precision.

My only question is the actual tolerance itself.  The value

eps * real (pivot)

works.  But, is there a theoretical justification for that number given the ichol algorithm?  I don't know because I don't know this math well enough to say.

If we don't have a theoretical basis, another way would be to try and determine the tolerance empirically.  Have the system print the value of

std::fabs(pivot.imag ())
std::fabs(pivot.real ())

and then try it on a few failing matrices.  We then set the comparison threshold to something that works for these matrices and see if there are further bug reports.

Rik <rik5>
Group administrator
Fri 10 Feb 2023 09:10:46 PM UTC, original submission:  

With some high (but "safe") optimization
(-O2 -mavx -mavx2 -flto) while using gcc 11 and gcc12
"test ichol" fails with

octave:1> test ichol
***** test
 opts.type = "nofill";
 opts.michol = "off";
 assert (nnz (tril (A1)), nnz (ichol (A1, opts)));
 assert (nnz (tril (A2)), nnz (ichol (A2, opts)));
 assert (nnz (tril (A3)), nnz (ichol (A3, opts)));
 assert (nnz (tril (A4)), nnz (ichol (A4, opts)));
 assert (nnz (tril (A5)), nnz (ichol (A5, opts)));

!!!!! test failed
ichol: non-real pivot encountered.  The matrix must be Hermitian.

In particular:

octave:2> A5 = [ 0.37, -0.05,         -0.05,  -0.07;
       -0.05,  0.116,         0.0,   -0.05 + 0.05i;
       -0.05,  0.0,           0.116, -0.05;
       -0.07, -0.05 - 0.05i, -0.05,   0.202];
 A5 = sparse (A5);
octave:4> ishermitian (A5)
ans = 1
octave:5> ichol(A5)
error: ichol: non-real pivot encountered.  The matrix must be Hermitian.
error: called from
    ichol at line 246 column 9

I think this is due to line 60 in _ichol_.cc

  if (pivot.imag () != 0)

If I replace it with

if (std::fabs(pivot.imag ()) > std::numeric_limits<double>::epsilon() * std::fabs(pivot.real ()))

test passes.

This applies both to stable and dev versions.


Dmitri A. Sergatskov <dasergatskov>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #54938:  bug63785.cset added by rik5 (1KiB - application/octet-stream)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by dasergatskov (Submitted the item)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.


    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-07-19 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
        Fixed ReleaseNone 8.3.0
    2023-07-18 rik5 Attached File- Added bug63785.cset, #54938
        StatusNone Patch Submitted

    Back to the top

    Powered by Savane 3.13-ee12.
    Corresponding source code