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:  None Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * 8.0.90 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment (Rich Markup)
   

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.
--

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.
--

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.
--

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

```cpp
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>
Project 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.
--






Dmitri A. Sergatskov <dasergatskov>

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

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 project members can vote.

     

     

     

     

    No changes have been made to this item

    Back to the top

    Powered by Savane 3.11