1 //------------------------------------------------------------------------------
2 // GB_mx_xsame64: check if two FP64 arrays are equal (ignoring zombies)
3 //------------------------------------------------------------------------------
4 
5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
6 // SPDX-License-Identifier: Apache-2.0
7 
8 //------------------------------------------------------------------------------
9 
10 #include "GB_mex.h"
11 
GB_mx_xsame64(double * X,double * Y,int8_t * Xb,int64_t len,int64_t * I,double eps)12 bool GB_mx_xsame64  // true if arrays X and Y are the same (ignoring zombies)
13 (
14     double *X,
15     double *Y,
16     int8_t *Xb,     // bitmap of X and Y (NULL if no bitmap)
17     int64_t len,    // length of X and Y
18     int64_t *I,     // row indices (for zombies), same length as X and Y
19     double eps      // error tolerance allowed (eps > 0)
20 )
21 {
22     if (X == Y) return (true) ;
23     if (X == NULL) return (false) ;
24     if (Y == NULL) return (false) ;
25     for (int64_t i = 0 ; i < len ; i++)
26     {
27         if (Xb != NULL && Xb [i] == 0)
28         {
29             // ignore X [i] and Y [i] if they are not in the bitmap
30             continue ;
31         }
32         // check X [i] and Y [i], but ignore zombies
33         if (I == NULL || I [i] >= 0)
34         {
35             int c = fpclassify (X [i]) ;
36             if (c != fpclassify (Y [i])) return (false) ;
37             if (c == FP_ZERO)
38             {
39                 // both are zero, which is OK
40             }
41             else if (c == FP_INFINITE)
42             {
43                 // + or -infinity
44                 if (X [i] != Y [i]) return (false) ;
45             }
46             else if (c != FP_NAN)
47             {
48                 // both are normal or subnormal, and nonzero
49                 double err = fabs (X [i] - Y [i]) / fabs (X [i]) ;
50                 if (err > eps) return (false) ;
51             }
52         }
53     }
54     return (true) ;
55 }
56 
57