Line data    Source code

       1             : """
       2             :     isapprox(A, B)
       3             :     A ≈ B
       4             : 
       5             : `isapprox(A, B)` for *sparse* matrices `A` and `B`. In contrast to the
       6             : standard implementation on `AbstractArray`s zeros at `A[i,j] == B[i,j]
       7             : == 0` are not taken into account. The algorithm is `O(nnz(A))`.
       8             : 
       9             : !!! note
      10             :     This method avoids testing for zero as in `norm(A - B) ≈ 0`, which
      11             :     is what the standard implementation does.
      12             : 
      13             : See also [`isapproxsymmetric`](@ref)
      14             : """
      15           4 : function Base.isapprox(A::SparseOp, B::SparseOp; kwargs...)
      16           2 :     size(A) == size(B) || return false
      17             : 
      18             :     # Note: Using nnz(spfilter(approxnonzero, A-B)) is tempting but requires
      19             :     #       x ≈ 0, which is difficult/imprecise. It is hard to get good defaults
      20             :     #       for atol (or atol from STOL).
      21             :     #       The implementation below is less efficient but more "accurate".
      22             : 
      23           2 :     i, j, v = spfind(A)
      24             : 
      25             :     # Inefficient but simple.
      26             :     # (Could provide a more efficient scan over coordinate forms.)
      27           2 :     for k in eachindex(i)
      28          60 :         isdifferent(k) =
      29             :             (!isapprox(v[k], B[i[k], j[k]]; kwargs...) &&
      30             :             !(abs(v[k]) < 10eps(typeof(v[k]))))
      31             : 
      32          20 :         isdifferent(k) && @show ((i[k], j[k]), v[k], B[i[k], j[k]])
      33          20 :         isdifferent(k) && return false
      34          20 :     end
      35           0 :     true
      36             : end
      37             : 
      38             : """
      39             :     isapproxsymmetric(A[; atol, rtol])
      40             :     A ≈ A'
      41             : 
      42             : Check if the sparse matrix `A` is approximately symmetric.
      43             : 
      44             : A typical use case is checking symmetry in the presence of numerical
      45             : inaccuracies due to different order of operations for entries in the
      46             : upper/lower triangles, e.g., `L = G'*A*G`.
      47             : 
      48             : The algorithm is `O(nnz(A))`.
      49             : 
      50             : See also [`isappox`](@ref)
      51             : """
      52          12 : function isapproxsymmetric(A::SparseOp;
      53             :                            atol::Real=0, rtol::Real=atol>0 ? 0 : sqrt(eps()))
      54           6 :     m, n = size(A)
      55             : 
      56           6 :     (m != n) && return false
      57           6 :     issymmetric(A) && return true
      58             : 
      59           6 :     _, _, err = spfind(A - A')
      60           6 :     _, _, v = spfind(A)
      61             : 
      62           6 :     tol = max(atol, rtol * norm(v) / n)
      63             : 
      64           6 :     norm(err) <= tol
      65             : end